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FOREWORD 


This  user's  manual  includes  new  features  of  the  STAGS  Computer  Code 
developed  by  the  Lockheed  Palo  Alto  Research  Laboratory,  Palo  Alto, 
California  under  sponsorship  of  the  Lockheed  Independent  Research 
Program,  NASA-Langley  Research  Center  (NASA/LRC)  (Contract  No. 
NAS1-10843)  and  the  Air  Force  Flight  Dynamics  Laboratory  (AFFDL) 
Contract  No.  F33615-69-C-1523).  The  AFFDL  Contract  was  admin¬ 
istered  under  the  Structures  Division,  with  Mr.  T.  N.  Bernstein 
(AFFDL/FBR)  as  Project  Engineer. 


This  report  was  completed  in  October  1972  and  covers  work  performed 
between  April  1969  and  October  1972.  The  supervision  of  this  project 
was  provided  by  Mr.  B.  O.  Almroth  of  the  Structural  Mechanics 
Laboratory,  LMSC. 


This  technical  report  has  been  reviewed  and  is  approved. 


Air  Force  Flight  Dynamics  Laboratory 
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ABSTRACT 

This  user's  manual  presents  STAGS,  a  comprehensive  computer  code. 
STAGS  is  intended  for  the  static  analysis  of  arbitrary  shells  Including 
the  effects  of  nonlinearities  caused  by  material  behavior  and  finite 
deformations.  Collapse  loads  based  on  nonlinear  analysis  ccn  be  com¬ 
puted  as  well  as  buckling  loads  based  on  classical  bifurcation  buckling 
theory  with  linear  prestress.  Arbitrary  thermal  and  mechanical  load¬ 
ings  can  be  specified.  The  manual  provides  instructions  for  use  of  the 
code  and  presents  sample  problems  and  solutions.  The  theoretical 
basis  for  the  program  also  is  presented.  The  program  is  under  devel¬ 
opment  and  the  version  presented  here  is  expected  to  be  updated  in  1973. 
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INTRODUCTION 


STAGS  is  a  computer  c«x1e  developed  to  analyze  the  behavior  of  general  shells  under 
arbitrary  static  thermal  and  mechanical  loading.  Nonlinearities  caused  by  material 
liuhavior  and  finite  deformations  are  accounted  for.  The  STAGS  analysis  is  l>used  on 
an  energy  formulation.  Derivatives  which  appear  in  the  energy  expression  are  replaced 
by  their  two-dimensional  finite  difference  approximations.  When  the  energy  is  rendered 
stationary,  the  result  is  a  system  of  nonlinear  algebraic  equations,  which  are  solved 
by  use  of  a  modified  Newton- Raphson  method. 

STAGS  is  an  outgrowth  of  work  on  the  buckling  of  cylindrical  panels  with  nonuniform 
membrane  stresses.  This  work  was  initiated  at  LMSC  in  1963  under  the  sponsorship 
of  NASA  Marshall  Space  Flight  Center  (Ref.  1).  The  basic  nonlinear  computer  program 
for  cylindrical  shells  with  cutouts  (Ref.  2)  and  a  linear  version  including  analysis  of 
free  vibrations  (Ref.  3)  were  developed  under  the  LMSC  Independent  Research  Program. 
Under  contract  with  the  Naval  Ship  Research  and  Development  Center  (NSRDC),  the 
linear  version  of  the  code  was  developed  to  include  shells  of  revolution  with  smooth 
but  otherwise  arbitrary  cutouts  (Ref.  4).  This  user's  manual  describes  the  latest  ver¬ 
sion  of  STAGS  which  includes  features  developed  during  the  last  three  years  in  separate 
studies  performed  under  the  sponsorship  of  the  Air  Force  Flight  Dynamics  Laboratory 

f 

(AFFDL),  the  Air  Force  Space  and  Missile  Systems  Organization  (SAMSO),  and  the 
LMSC  Independent  Research  Program.  In  the  AFFDL  effort,  the  nonlinear  capability 
was  extended  to  shells  of  more  general  shape  and  with  cutouts  of  arbitrary  contour. 

In  addition,  inelastic  deformations  were  Introduced  and  a  capability  to  handle  a  finite 
difference  grid  with  variable  nodal  point  spacing  was  added.  In  a  parallel  effort  spon- 
sored  under  the  LMSC  Independent  Research  Program,  the  equations  were  generalized 

. 

to  include  nonorthogonal  coordinates  (Ref.  5).  Further  expansion  of  STAGS  was  ac- 

[ 

complished  under  the  SAMSO  sponsored  study.  Provisions  were  made  to  allow  both 

f 
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temperature  and  material  properties  to  vary  over  the  surface  and  through  the  thickness 
of  the  shell.  A  bifurcation  buckling  branch  was  added  and  parameter  studies  were  made 
to  evaluate  the  applicability  of  bifurcation  buckling  analysis  for  shells  of  general  shu|>e. 


This  user's  manual  described  features  developed  under  each  of  these  studies.  For 
further  information  concerning  the  analyses  and  parameter  studies,  the  reader  is 
referred  to  final  reports  issued  under  these  studies  (Refs.  6  and  7). 

This  report  begins  with  a  discussion  of  collapse  and  bifurcation  buckling.  Then,  a 
description  of  the  scope  and  limitations  of  the  code  is  given  followed  by  details  of  the 
analysis.  Next,  the  information  necessary  to  use  STAGS  is  presented. 


x 


Section  1 

COLLAPSE  OF  SHELL  STRUCTURES 


A  stress  analyst  usually  defines  the  buckling  load  as  the  value  of  the  load  at  which  the 
fundamental  branch  of  the  load -displacement  curve  is  intersected  by  a  branch  corre¬ 
sponding  to  a  buckled  equilibrium  form.  This  definition  is  meaningful  if  the  shell 
possesses  a  high  degree  of  symmetry  as  in  the  case  of  axisymmetrically  loaded  shells 
of  revolution.  In  reality,  small  deviations  from  the  nominal  geometry  will  destroy 
this  symmetry.  In  this  case,  bifurcation  in  the  equilibrium  curve  does  not  occur,  and 
the  classical  concept  of  buckling  exists  only  as  an  idealization. 

It  also  may  be  noticed  that  the  behavior  exhibited  by  the  shell  as  the  buckling  load  is 
reached  varies  drastically  from  one  case  to  another.  As  an  example,  buckling  of  a 
spherical  cap  under  uniform  pressure  corresponds  to  a  complete  collapse  of  the  shell. 
However,  if  the  same  spherical  cap  is  subjected  tb  a  radial  point  force,  buckling  de¬ 
termined  through  a  bifurcation  analysis  may  correspond  to  a  load  level  at  which  a 
gradual  change  in  deformation  pattern  starts  to  develop.  Loss  of  stability  of  the  funda¬ 
mental  equilibrium  configuration  then  corresponds  only  t_  an  imperceptible  change  in 
the  stiffness  of  the  shell  and  is  of  no  consequence  to  the  designer. 

In  these  cases,  a  nonlinear  analysis  of  the  axisymmetric  prebuckling  configuration  is 
combined  with  a  buckling  analysis  based  on  the  assumption  that  the  nonsymmetric  dis¬ 
placements  are  infinitesimal.  A  two-dimensional  nonlinear  analysis  gives  a  more 
complete  description  of  the  shell  behavior  but  is  much  more  expensive.  Hence,  for 
shells  with  a  high  degree  of  symmetry,  the  classical  type  of  buckling  analysis  will 
remain  an  important  instrument.  Buckling  in  the  classical  sense  means  that  the 
equilibrium  of  the  fundamental  branch  becomes  unstable  relative  to  some  deformation 
mode  which  is  orthogonal  to  the  prebuckling  pattern.  For  a  shell  of  a  general  shape, 
it  can  be  expected  that  all  deformation  modes  are  present  in  the  displacement  pattern 
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corresponding  to  the  fundamental  branch.  In  this  case,  a  complete  nonlinear  analysis 
is  the  only  rigorous  solution  to  the  problem. 

If  a  plane  of  symmetry  exists  with  respect  to  loading  as  well  as  geometry,  bifurcation 
is  possible  into  modes  which  are  antisymmetric  relative  to  this  plane.  Thus,  a  non¬ 
linear  analysis  which  uses  symmetry  to  reduce  the  amount  of  calculations  should  be 
supplemented  by  a  bifurcation  analysis  with  respect  to  antisymmetric  displacements. 
Also,  it  is  quite  possible  that  a  classical  buckling  analysis  may  give  good  estimates 
fo*  the  collapse  load  even  in  cases  for  which  it  is  not  strictly  applicable.  This  would 
be  the  case  if  little  change  in  load  distribution  or  shell  geometry  occurs  before  a  fatal 
type  of  buckling  pattern  starts  to  develop.  For  example,  if  a  short  cylinder  is  loaded 
by  nonuniform  external  pressure,  the  bifurcation  analysis  may  give  a  good  approxima¬ 
tion  of  the  collapse  load,  but  it  may  not  represent  an  adequate  solution  for  a  longer 
cylinder. 
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Section  2 

SCOPE,  LIMITATIONS,  PITFALLS 


The  computer  program  STAGS  (Structural  Analysis  of  General  Shells)  performs  a  non¬ 
linear  analysis  of  shells  by  use  of  a  two-dimensional  finite  difference  approach.  Dis¬ 
placement  and  stress  histories  are  computed  corresponding  to  a  given  history  of 
applied  load,  displacement,  or  temperature.  Two  branches  of  STAGS  ar^  described 
in  this  manual.  The  first  branch  is  restricted  to  elastic  material  behavior  and  in¬ 
cludes  three  sub-branches:  (1)  nonlinear  collapse  analysis,  (2)  linear  analysis,  and 
(3)  buckling  analysis  based  on  the  classical  bifurcation  approach  with  a  linear  pre¬ 
buckling  analysis.  Collapse  loads  are  fount'  as  limit  points  in  the  nonlinear  load 
displacement  curve.  The  program  also  is  useful  for  postbuckling  analysis  of  shells 
which  behave  according  to  classical  buckling  theory  and  for  studies  of  the  influence 
of  Imperfections. 

The  second  branch  is  for  inelastic  material  behavior  but  does  not  allow  temperature 
or  material  properties  to  vary  with  any  of  the  space  coordinates.  Only  the  nonlinear 
analysis  sub-branch  is  included  here. 

STAGS  applies  to  any  shell  for  which  a  reference  surface  and  a  suitable  set  of  gridlines 
can  be  mathematically  defined.  In  general,  the  user  of  the  program  provides  a  sub¬ 
routine  describing  the  geometry,  but  several  such  routines  for  standard  geometries 
are  permanently  Included  in  the  program  (these  are  listed  in  Section  6).  For  the 
elastic  branch  the  shell  wall  thickness  can  be  varied,  and  elastic  properties  are 
allowed  to  vary  with  the  shell  coordinates  and  through  the  thickness.  Cutouts  in  the 
shell  wall  and  discrete  eccentric  stiffeners  are  included.  The  program  is  also  general 
relative  to  boundary  conditions  and  to  loading.  The  loading  can  be  applied  in  terms  of 
variable  surface  tractions,  point  forces,  or  lineloads.  Displacements,  such  as  uniform 
end  shortening  of  a  cylindrical  shell,  can  be  applied  if  desired  rather  than  fixed  loading, 
and  provision  is  made  for  thermal  loading. 
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The  general  equations  on  which  the  analysis  is  based  are  given  in  Section  3. 1  "Basic 
Equations,  "  3.  3  "Bifurcation,  "  and  3.  4  "Plasticity.  "  The  finite  difference  expressions, 
valid  for  a  variable  grid  spacing,  are  discussed  in  Section  3.5.  The  introduction  of  the 
finite  difference  approximations  into  the  energy  expression  leads  to  a  nonlinear  algebraic 
eq  !  tion  system.  Section  3.2  discusses  the  method  of  solution  of  this  system. 

Stiffeners  and  cutout  edges  must  follow  coordinate  lines,  or  rather,  the  coordinate  lines 
must  be  chosen  so  that  they  follow  boundaries,  internal  or  external,  and  the  direction  of 
stiffeners.  This  is  not  a  severe  program  limitation  because  the  capability  of  handling 
nonorthogonal  grids  has  been  included.  However,  analysis  of  more  complicated  shells 
will  require  some  user  skill. 

Torsional  stiffness  and  the  resistance  against  rotation  are  included  in  the  3tifiener 

strain  energy,  but  the  energy  caused  by  change  of  curvature  in  the  plane  of  the  shell 

has  been  omitted.  Note  that  shell  stiffeners  usually  are  of  the  type  shown  in  Fig.  2-1; 

l,  e.  ,  they  are  designed  to  minimize  bending  of  the  shell  surface.  Because  of  the 

bending  flexibility  of  the  thin  web,  the  line  of  intersection  between  the  frame  web  and 

the  shell  surface  can  twist  and  bend  (in  the  shell  plahe)  without  much  resistance  from 

the  stiffener.  In  such  a  case,  t  is  recommended  that  corresponding  stiffnesses  (J  and 

I  )  are  read  in  as  zero.  Then,  in  the  stress  output,  the  very  small  stresses  caused 
z 

by  the  changes  of  curvature  in  the  p^e  of  the  shell  will  be  ignored.  For  accurate 

analysis  of  stiffeners  whose  cross- 
section  is  likely  to  distort,  it  would  be 
necessary  to  consider  such  stiffeners  as 
shell  branches. 

The  input  format  requires  that  four  shell 
edges  be  specified.  For  a  shell  that  is 
closed  in  one  direction,  such  as  a  com¬ 
plete  shell  of  revolution,  the  program 
will  assume  that  the  shell  is  closed  if 
appropriate  boundary  conditions  are 
specified  on  opposite  edges.  For  all 
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the  special  geometry  routines  presently  in  the  program,  the  shells  may  be  closed 
only  in  the  tj  direction. 

On  each  of  the  boundaries,  the  input  parameters  can  be  used  to  specify  restraint  on 
either  of  the  three  displacement  components  or  on  the  rotation  about  the  tangent  to  the 
edge.  If  displacement  restraint  is  not  specified  for  one  or  more  of  these  quantities, 
the  analysis  will  be  based  on  appropriate  natural  boundary  conditions  (stress  free). 

If  more  complicated  boundary  conditions  are  used,  it  will  be  necessary  to  modify  the 
program  in  an  area  which  is  not  easily  accessible  to  the  user.  However,  displace¬ 
ment  restraints  off  the  bourd-ry  lines  can  be  introduced  by  using  specified  displace¬ 
ments  in  the  load  routine. 

The  most  severe  limitation  of  the  program  in  practical  analysis  is  dictated  by  com¬ 
puter  economy.  A  large  number  of  unknowns  is  generally  unavoidable  in  a  two  dimen¬ 
sional  numerical  analysis  regardless  of  the  solution  method  used.  Because  of  the 
coupling  between  the  unknowns  in  two  directions,  the  bandwidth  of  the  resulting  equation 
system  is  necessarily  large  in  comparison  to  the  bandwidth  for  one-dimensional  prob¬ 
lems.  In  addition,  the  nonlinearity  of  the  formulation  requires  an  iterative  procedure 
in  which  the  basic  linear  equation  system  is  solved  a  number  of  times.  Particularly  if 
plasticity  is  included,  the  load  step  size  must  be  small,  and  the  convergence  may  be 
rather  slow. 

Two  items  are  of  special  importance  relative  to  the  problem  of  computer  runtime. 

One  is  the  transformation  of  the  structure  into  a  model  which  is  suitable  for  analysis; 
the  other  is  the  choice  of  strategy  in  the  nonlinear  analysis.  Often,  it  is  impossible 
to  model  the  structure  in  a  straightforward  manner  and  considerable  engineering  skill 
may  be  needed.  Eventually  a  few  different  models  must  be  analyzed,  all  reflecting 
different  types  of  local  behavior.  The  strategy  in  the  analysis  involves  choice  of 
such  items  as  stepsize  and  convergence  criterion,  and  also  involves  the  use  of  initial 
imperfections  and  the  use  of  the  results  from  the  application  of  the  bifurcation  anal¬ 
ysis.  The  choice  of  a  proper  strategy  is  very  important  for  computer  economy;  this 
is  discussed  in  detail  in  Section  4. 
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For  maximum  generality  of  the  computer  program  it  is  necessary  to  introduce  basic 
data  (loads,  temperature,  thickness)  through  user-written  subroutines.  Such  routines 
are  discussed  in  Section  5.  Usually,  they  are  very  easy  to  derive,  but  geometry 
routines  may  be  complicated.  Because  of  the  importance  of  the  strategy  and  of  the 
modeling  of  the  structure,  it  is  concluded  that  considerable  user  skill  is  required  to 
achieve  the  full  potential  of  the  STAGS  program. 

The  basic  problem  dimensions  (number  of  nodal  points  in  each  direction)  are  restricted 
only  by  the  availability  of  mass  storage.  Some  program  limitations  are  imposed  by 
dimension  statements.  In  each  direction,  there  can  be  as  many  as  80  stiffeners  but 
only  of  30  different  types.  The  number  of  points  for  integration  through  the  thickness 
must  be  an  odd  number  and  it  may  not  exceed  9.  In  the  plasticity  analysis,  as  many 
as  10  material  components  may  be  used. 

The  unwary  user  may  encounter  pitfalls.  For  example,  use  of  a  too  coarse  grid  may 
lead  to  inaccurate  results.  This  problem  is  discussed  in  Section  4.  Another  problem 
arises  when  a  buckling  mode  starts  to  develop,  and  its  amplitude  is  small  in  compari¬ 
son  to  the  total  displacement.  The  convergence  criterion  then  may  not  be  sufficiently 
sharp  to  catch  the  growth  of  this  new  deformation  mode.  This  problem  is  discussed  in 
Section  4.  In  particular,  it  is  possible  that  the  shell  may  buckle  in  a  mode  which  is 
antisymmetric  with  respect  to  a  plane  about  which  symmetric  behavior  has  been  assumed. 
Such  an  occurrence  probably  will  be  revealed  through  an  inspection  of  the  buckling  mode. 
In  the  bifurcation  analysis,  it  is  possible  to  use  different  boundary  conditions  for  pre¬ 
buckling  and  incremental  displacements. 

There  is  also  a  possibility  that  the  problem  as  defined  by  the  program  user  is  not  well 
posed.  Difficulties  will  arise  if  the  boundary  conditions  allow  rigid  body  displacements. 
This  will  make  the  system  extremely  ill  conditioned;  it  would  be  singular  except  for  the 
truncation  errors  in  the  finite  difference  expressions.  The  difficulty  will  occur  only  in 
analysis  of  a  free  body  which  is  subjected  to  a  self-equilibrating  force  system.  Rigid 
body  displacements  can  be  restrained  by  specifying  displacement  constraints  as  discussed 
in  Section  6  (input  description,  L-l  and  L-2  cards).  Because  of  the  aforementioned 
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truncation  error,  a  small  force  may  develop  at  restrained  points.  Therefore,  it  is 
advisable  to  fix  points  at  which  the  structure  is  relatively  stiff.  Results  would  be 
meaningless  if  a  load  is  defined  on  an  edge  at  which  the  displacements  are  restrained 
in  the  direction  of  the  load.  The  same  would  be  the  case  if  a  prescribed  displacement 
is  in  conflict  with  displacement  boundary  conditions  on  an  adjacent  edge. 

For  a  general  shell  shape,  the  bifurcation  buckling  branch  of  the  program  does  not  cor¬ 
respond  to  a  rigorous  application  of  stability  theory.  The  results  may  or  may  not 
represent  good  approximations.  This  problem  is  discussed  in  Ref.  3. 
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Section  3 
ANALYSIS 


3. 1  BASIC  EQUATIONS 

For  shells  of  a  mote  general  shape  than  the  shell  of  revolution,  it  is  not  possible  to 
separate  the  governing  differential  equations.  The  analysis  thus  requires  the  use  of 
two  Independent  space  variables,  and  the  numerical  analysis  is  drastically  encumbered. 
In  addition,  if  the  collapse  load  of  the  shell  is  to  be  determined,  a  nonlinear  analysis 
is  encountered.  The  STAGS  program  is  based  on  a  discretization  of  the  total  potential 
energy  by  use  of  finite  difference  expressions. 

The  energy  method  was  used  because  it  simplifies  the  handling  of  shell  cutouts  and 
discrete  stiffening.  For  a  general  shell,  the  surface  coordinates  are  chosen  for 
practical  reasons  to  coincide  with  shell  boundaries.  In  this  case,  they  are  not  neces¬ 
sarily  lines  of  curvature,  and  the  basic  equations  must  be  written  in  terms  of  non- 
orthogonal  coordinates.  For  brevity,  the  basic  equations  are  given  here  in  tensorial 
form,  but  a  formulation  in  terms  of  physical  components  is  available  in  Bef.  2. 

The  expression  for  the  strain  energy  is 


1  E 

U  =  2 - 2 

2  1  -  v2 


( 1  -  r. )  ft  *  fj  %,  j 


<D 


The  expressions  for  strains  and  curvature  changes  are 


eap  ~  2^ap  +  +  2  ^a^p  +  2  ^pa^.P 

'af-tya*  *6  V  -  bc^% 


(2) 
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where  and  pa  are  the  displacement  gradients  defined  by 


^ap  ua  I  p  ^ap  w 

f>a  '  »•  a  +  >4  ap  <3> 

and  the  curvature  tensor  b  „  and  the  normal  displacement  w  are  defined  with  respect 
to  the  inner  shell  normal  vector.  The  curvature  tensor  differs  from  that  by  Sanders 
(Ref.  8)  because  it  is  valid  for  larger  out-of -plane  rotations.  A  complete  derivation  is 
given  in  Ref.  9. 

3.2  SOLUTION  METHOD 

The  numerical  solution  is  based  upon  a  two-dimensional  finite  difference  approximation. 
The  shell  surface  is  covered  with  mesh  lines  parallel  to  the  coordinate  lines,  and  the 
freedoms  of  the  system  are  the  normal  displacements,  w ,  at  the  grid  points  and  the 
tangential  displacements,  u  and  v,  at  points  between  adjacent  gridpoints. 

After  replacement  of  the  displacement  functions  and  their  derivatives  in  the  governing 
equations  by  finite  difference  approximations  (see  Section  3.  5),  the  strain  energy  den¬ 
sity  at  mesh  station  i  can  be  written  in  the  form 

AU1  *=  |  Z1*  Dl  Z*  (4) 

where  D*  is  a  6  x  6  positive  definite  matrix  of  constants  and  Z*  is  a  column  vector 
of  strain  and  curvature  changes  at  station  i .  D*  and  Z*  are  functions  of  the  geo¬ 
metric  parameters  of  the  shell;  in  addition,  D*  is  dependent  on  the  material  properties. 
Z1  is  a  nonlinear  (quadratic)  function  of  the  displacement  unknowns  and  thus  AU*  is 
a  fourth-order  polynomial.  The  vector  of  stress  resultants  S1  at  station  i  is  given 
by 


S*  =  D1  Z* 


(5) 
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The  total  potential  energy,  V,  of  the  shell  is  obtained  by  combination  of  the  strain 
energy  and  the  work  done  by  the  external  forces. 

V  =  U  -  W  /( 


where 


m  .  . 

U  ■  £  Aljl  '  a* 
i 


and 


W  =  X*  •  F 

Here  X  denotes  the  vector  of  displacement  components,  F  is  the  vector  of  external 
forces  and  a .  is  the  area  of  the  1th  subregion.  A  necessary  condition  for  static 
equilibrium  is  that  the  total  potential  energy  be  stationary.  This  condition  requires 
the  vanishing  of  the  first  variation  of  V  and  leads  to  the  equation 

LX  *=  F  (7) 


where  the  operator  L  is  defined  by 

LX  =  Grad  U  (8) 

L  is  thus  a  "stiffness"  operator  which  relates  displacement  components  and  external 
forces  and  is  nonlinear  in  the  general  case. 

When  only  linear  terms  are  include^  in  the  definition  of  the  strains  and  changes  in 
curvature,  L  is  a  linear  operator  which  may  be  readily  represented  in  matrix  form 

[see  Section  3.3  and  Eqs.  (22)  and  (23)] .  In  this  case,  the  matrix  is  positive  definite 
(with  the  choice  of  proper  boundary  conditions)  and  Eq.  (7) 
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LX 


may  be  solved  by  one  of  many  direct  or  iterative  methods.  However,  when  geometric 
nonlinearities  (i.e. ,  rotations)  are  included,  L  becomes  a  polynomial  operator  of 
third  degree  and  iterative  methods  must  be  employed  for  solution  of  the  equations. 

For  a  general  collapse  analysis,  it  is  necessary  to  solve  the  operator  equation, 

Eq.  (?)  for  a  sequence  of  applied  loads.  In  fact,  the  only  practical  method  consists 
of  a  sequence  of  load  steps  chosen  so  that  the  initial  solution  is  nearly  linear  and  sub¬ 
sequent  solutions  change  only  moderately  from  one  step  to  the  next.  Such  a  procedure 
is  mandatory  for  two  reasons:  first,  the  feasibility  of  the  iterative  methods  of  solu¬ 
tion  depend  on  reasonably  good  initial  approximations  and  second,  a  reliable  detection 
of  collapse  requires  such  a  stepwise  procedure  because  of  the  non-uniqueness  of 
solutions  to  nonlinear  equation  systems. 

Thus,  at  the  1th  load  step  in  a  collapse  analysis,  the  operator  equation 


LX 


(6) 


must  be  solved  where  is  the  vector  of  constants  generated  by  the  applied  load 
(mechanical,  thermal,  plastic  pseudo  loads). 


A  brief  description  of  Newton's  method  and  the  modified  Newton  method  for  the  case 
of  a  function  of  one  variable  displays  the  principal  features  of  the  methods  used  for 
the  solution  of  the  nonlinear  equation  in  STAGS. 

Newton's  method  for  the  solution  of  the  problem  g(x)  =  0  (see  Fig.  3-1)  is  defined  by 


Vl  ■  xn  - 


(») 


The  iteration  converges  quadratically  to  the  solution  provided  that  the  initial  estimate 
Xq  is  sufficiently  accurate.  A  geometric  Interpretation  of  Newton's  method  is  well 
known;  from  any  point  on  the  curve  g(x),  the  next  approximation  is  obtained  by  extend¬ 
ing  the  tangent  to  the  curve  at  the  point  to  the  x-axis.  In  Fig.  3-2,  the  "modified"  Newton 
iteration  is  Illustrated.  The  method  is  defined  by  the  equation 
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Fig.  3-1  Newton’s  Method 


Fig.  3-2  Modified  Newton  Method 
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xn+  1  =  x„  - 


(10) 


Geometrically,  this  iteration  corresponds  to  extending  lines  from  the  curve  g(x)  to 
the  axis  which  are  all  parallel  to  the  tangent  at  xq  .  While  the  convergence  of  the 
modified  Newton  method  obviously  is  slower  than  that  of  the  standard  Newton  method, 
it  avoids  repeated  computation  of  g'(x).  The  most  effective  strategy  usually  calls 
for  periodic  recomputation  of  g'(x). 

Both  Newton's  method  and  the  modified  Newton  method  have  been  generalized  to  n- 
dimensional  Euclidean  space. 


The  treatment  of  the  complete  nonlinear  solution  of  Eq.  (7)  as  well  as  of  bifurcation 
buckling  is  facilitated  by  introduction  of  the  concept  of  the  derivative  L'  of  L  (Ref.  6). 
In  particular,  for  the  operator  L,  the  derivative  L'  (sometimes  called  the  Frechet 
derivative  of  L )  is  an  n  by  n  matrix  whose  elements  are 


a2u 


9x(i)ax(J) 


(11) 


Most  of  the  elementary  properties  of  ordinary  derivatives  also  hold  for  the  Frechet 
derivatives  L'  of  an  operator  L. 


Because  L'  is  a  function  of  a  particular  displacement  vector  X  (unless  the  nonlinear 
terms  are  dropped)  the  Frechet  derivative  will  usually  be  denoted  to  indicate  this 
dependence.  With  the  use  of  the  derivative  L'  of  the  operator  L,  Newton's  method 
may  be  readily  generalized  to  obtain  a  solution  of  Eq.  (7).  The  iteration  is  defined  by 


IVxk*i-v  =  F-“k  < 

If  Xq  is  sufficiently  close  to  a  solution  X  and  if  is  not  a  singular  matrix,  the 
iteration  converges  to  X .  Under  these  assumptions,  it  also  can  be  shown  that  the 
converged  solution  is  unique  in  some  neighborhood  of  X  (Ref.  10). 


3-6 


Similarly,  with  the  aid  of  the  derivative  L^, ,  the  modified  Newton  method  may  be 
applied  to  the  operator  equation  Eq.  (  8  ).  The  general  form  of  the  iteration  then 
becomes 


He  <Xn+l'Xn>  =  Fi  '  <13> 

m 

This  iteration  formula  contains,  as  special  cases,  the  most  commonly  used  iterative 
methods.  For  example,  if  Xm  =0,  is  identical  with  the  matrix  A  corre¬ 

sponding  to  the  linearization  of  the  operator  L.  If  B  denotes  that  part  of  the  operator 
L  obtained  by  variation  of  the  3rd  and  4th  degree  terms  of  V ,  the  iteration  may  be 
written 


AX  .  *=  F.  -  BX  (14) 

n+1  in  '  ' 

This  Iteration  corresponds  to  repeated  solution  of  the  linear  equation  system  which  is 
obtained  by  substitution  of  the  previous  Iterates  for  the  unknowns  in  the  nonlinear  terms. 
The  method  is  often  effective  in  cases  for  which  the 'effects  of  nonlinear  terms  are 
small.  Unfortunately  in  practical  cases  (e.g. ,  a  cylinder  with  a  cutout),  the  iteration 
often  falls  to  converge  well  before  the  collapse  load  is  reached.  Furthermore,  when 
the  method  is  applicable,  a  bifurcation  analysis  would  frequently  be  an  even  more 
effective  choice. 

If  only  one  iteration  of  Eq.  (13)  is  performed  for  each  load  step  and  if  m  =  i ,  the 
resulting  method  is  equivalent  to  the  frequently  used  "Incremental  analysis. "  Both 
this  method,  and  the  standard  Newton  method,  require  the  recomputation  of  the  deriva¬ 
tive  matrix  and  its  factorization  one  or  more  times  for  each  load  step.  For  the 
"incremental  analysis, "  the  load  steps  must  be  chosen  small  enough  to  ensure  accuracy, 
but  the  Iteration  does  not  itself  provide  the  information  necessary  for  such  a 
determination. 

In  contrast  to  the  iterative  methods  described  above,  the  modified  Newton  method 
provides  accurate  solutions  Independent  of  the  size  of  the  load  step  (numerical  error 


3-7 


does  not  accumulate)  and  at  the  same  time  avoids  the  necessity  of  frequent  recompu¬ 
tation  and  factorization  of  the  derivative  matrix  L' .  This  latter  feature  exploits  a 
fundamental  advantage  of  the  finite  difference  treatment  used  here  which  permits  ex¬ 
tremely  rapid  evaluation  of  LXn  .  Even  when  the  best  available  methods  of  factoriza¬ 
tion  are  employed,  the  time  required  for  a  single  evaluation  of  LXr  is  many  times 
less  than  the  factorization  time.  The  modified  Newton  method  is  employed  in  STAGS 
to  combine  rigorous  results  with  the  most  economical  computational  effort.  The 
effective  use  of  the  modified  Newton  method  requires  choices  both  as  to  the  size  of 
load  steps  and  as  to  when  the  derivative  L'  should  be  recomputed  and  refactored. 

The  STAGS  program  contains  as  much  built-in  decision  making  capability  regarding 
these  questions  as  appears  feasible.  However,  it  is  still  necessary  for  the  user  of 
the  program  to  consider  the  best  overall  "strategy"  relating  to  these  choices  (see 
Section  4). 

3.3  BIFURCATION 

Note  that  the  mathematical  characterization  of  bifurcation  buckling  also  is  provided  by 
the  generalized  Newton  method.  Let  Xq  be  a  solution  of  Eq.  (7)  under  a  given 
vector  F  of  external  forces.  If  every  neighborhood  of  Xq  contains  another  vector  Y 
which  satisfies  the  equation 


LY  =  F 


(1?) 


then  bifurcation  is  said  to  take  place  for  the  shell  under  the  load  F.  From  the  previous 
remarks  on  the  conditions  for  convergence  of  Newton's  method  to  a  unique  solution, 
it  follows  that  a  necessary  condition  for  bifurcation  is  that  LL  be  a  singular  matrix, 


det  (L*  )  =0  (16) 

0 

Classical  bifurcation  buckling  theory  may  be  obtained  easily  from  Eq.  (16).  It  is 
assumed  that  may  be  written 
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(17) 


X0  "  AXL 


where  is  the  linear  solution  for  a  load  vector  F^.  Thus,  Eq.  (16)  becomes 


det(^xL)  -  0 


(18) 


Equation  (18)  is  an  algebraic  eigenvalue  problem  of  the  form 


det  (A  -  XB  -  X  C)  =  0 


(19) 


In  classical  bifurcation  theory,  the  C  matrix,  which  arises  from  the  prebuckling  rota¬ 
tions,  is  often  omitted  and  the  eigenvalue  problem 

AX  «  XBX  (20) 


is  obtained. 

When  bifurcation  exists  but  the  prebuckling  displacements  are  not  linear,  the  solution 
of  Eq.  (16)  generally  requires  a  stepwise  procedure.  One  such  method  is  given  by  the 
equations 


det  (U  X.)  =  0  (21) 

Ak+  1  K 

*k+l  “*k+l*k 

where  Xq  is  the  linear  solution  A  sequence  of  eigenvalue  problems  is  solved  and, 
if  the  method  is  successful,  X^  approaches  one.  A  nonlinear  bifurcation  treatment, 
equivalent  to  Eq.  (21)  was  presented  in  Ref.  11  and  has  been  used  successfully  to  study 
a  large  variety  of  problems.  For  the  two-dimensional  problems  under  consideration 
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here,  it  appears  that  such  methods  may  be  nearly  as  costly  as  the  complete  nonlinear 
analysis  available  in  STAGS.  Consequently,  only  a  classical  bifurcation  buckling 
analysis  is  implemented  in  the  STAGS  program. 

The  formation  of  the  A  and  B  matrices  of  Eq.  (19)  will  be  considered  briefly.  The 
elements  of  the  Frechet  derivative  matrix  L^Xl  (whi°h  define  the  matrices  A  and  B) 
are  determined  according  to  Eq.  (11).  The  rules  for  computing  derivatives  of  poly¬ 
nomials  are  easily  programmed,  and  the  formation  of  the  A  and  B  matrices  therefore 

is  well  suited  to  automatic  treatment  on  the  computer.  Thus,  for  example,  if  X... 

th  th  ' 1 > 

and  X^  j  are  the  i  and  j  displacement  components,  by  use  of  Eq.  (4),  (5)  and 

(6)  the  following  is  obtained: 


92U 


9X(i)9X(l) 


m  2  k 

V  k  a  A  u 

Z  a  ax...  ax 

k=i  '■ 


The  k' 


th 


term  of  this  sum  is 


1*41^ 

3X(i)9X(i) 


g2  zk* 


+ 


k*  k 

az  nk  az 

9X(t)  9X(J) 


(22) 


(23) 


In  the  first  term  on  the  right  hand  side  of  Eq.  (23),  note  that  S**  is  the  linear  stress 

Lj 

resultant  vector  at  station  k  and  that  only  the  quadratic  terms  (rotations)  need  be  con¬ 
sidered  in  forming  the  partial  derivatives 


a2  Zk* 

9X<i)  9X(J) 


Contributions  from  this  term  go  into  the  B  matrix.  Assuming  the  prebuckling  rotations 
may  be  neglected  for  the  classical  theory,  the  last  term  of  Eq.  ( 23)  generates  additions 
only  to  the  A  matrix.  The  A  matrix  then  is  identical  with  the  linear  stiffness  matrix. 
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If  the  prebuckhng  rotations  are  included  to  treat  nonlinear  bifurcation,  the  last  term 
of  Eq.  (23)  generates  a  C  matrix  and  provides  additional  contributions  to  the  B  matrix. 
In  this  case,  the  prebuckling  stress  resultant  vector  S  would  also  include  nonlinear 
terms. 

In  conclusion,  it  should  be  noted  that  the  use  of  operation  notation  can  be  (and  often  is) 
avoided.  For  example,  a  nonlinear  equation  system  can  be  obtained  by  writing  the  dis¬ 
placements  at  step  n  +  1  as  the  sum  of  the  known  displacements  at  step  n  plus  an 
increment  of  displacement.  Newton’s  method  can  be  derived  for  functions  of  many 
variables  by  considering  Taylor  series  expansions  (see  Ref.  12).  Bifurcation  buckling 
theory  may  be  based  on  the  theory  of  adjacent  equilibrium  states.  In  all  of  these  cases, 
the  development  is  rather  complicated  and  burdened  with  excessive  detail.  The  use  of 
operator  methods,  however,  permits  the  immediate  application  of  a  well  developed 
theory.  The  concise  operator  notation  facilitates  manipulation  and  programming  for 
the  computer.  In  particular,  the  definition  and  computation  of  the  basic  matrices  of 
Eq.  (19)  are  greatly  simplified.  The  recipe  is  outlined  in  Eq.  (22)  and  (23)  and  can  be 
performed  using  straightforward  algebraic  procedures.  Finally,  the  relations  between 
a  complete  nonlinear  analysis ,  linear  (classical)  and  nonlinear  bifurcation  theory  and 
Newton's  method  are  clarified. 

3.4  PLASTICITY 

Introduction  of  inelastic  behavior  has  been  done  within  the  framework  of  the  energy 
principle  on  which  the  elastic  analysis  was  based.  Essentially,  the  plastic  deforma¬ 
tions  are  considered  as  load  terms;  they  are  completely  analogous  to  thermal  expan¬ 
sions  except  that  they  are  not  known  in  advance.  A  series  of  elastic  problems  are 
solved  by  use  of  the  energy  principles  in  which  the  "load  terms"  are  gradually  modi¬ 
fied  until  they  correspond  to  the  computed  state  of  stress  and  to  specified  nonlinear 
stress  strain  relations  at  all  points  over  the  shell  surface  and  through  the  shell 
thickness. 

The  plasticity  theory  used  has  been  proposed  by  Besseling  (Ref.  13),  and  is  based 
on  a  principle  which  originally  was  suggested  by  White.  This  theory  is  very  promising 
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because  it  is  rather  simple  in  its  application  yet  retains  such  features  as  strain  harden¬ 
ing  and  the  Bauschinger  effect. 

The  White- Besseling  theory,  as  applied  here,  assumes  that  the  material  consists  of 
several  components  which  have  identical  elastic  properties  and  exhibit  ideal  plasticity 
(no  strain  hardening)  but  have  different  yield  strength.  As  the  strain  is  the  same  in  all 
components,  the  stress-strain  curve  will  experience  a  decrease  in  slope  as  the  stress 
reaches  the  yield  limit  for  any  one  of  the  components;  the  respective  components  then 
cease  to  take  additional  load.  The  composite  thus  exhibits  strain  hardening  with  a 
piecewise  linear  stress-strain  relation.  Use  of  only  one  component  will,  of  course, 
result  in  application  of  ideal  plasticity  theory.  If  the  stress  is  reversed  after  loading 
beyond  the  yield  limit  for  one  or  more  components,  yield  will  occur  in  the  reversed 
direction  at  an  average  stress  in  the  composite  which  is  lower  than  the  stress  for 
original  yield.  This  is  demonstrated  in  the  uniaxial  stress -strain  curve  shown  in 
Fig.  3-3.  Tension  is  first  applied,  OAB,  beyond  the  yield  limit  and  followed  by  strain 
reversal,  BCD,  into  the  zone  of  yield  in  compression  .  The  yield  ellipse  for  the 
weakest  component  and  the  loading  history  in  this  component  are  also  shown  in  this 
figure.  Clearly,  yield  in  compression  will  occur  when  the  total  strain  is  -  2e^), 
i.e.  ,  the  yield  in  compression  occurs  at  a  much  lower  stress  if  the  material  pre¬ 
viously  has  been  subjected  to  tension  stresses  above  the  yield  point.  To  introduce  the 
Bauschinger  effect  this  way  is  appealing  because  it  reflects  the  microstress  theory 
which  now  generally  is  accepted  as  the  explanation  of  the  Bauschinger  effect. 

The  White -Besseling  plasticity  theory  is  implemented  in  the  computer  program  in  the 
following  manner: 

(1)  The  inelastic  behavior  of  the  material  is  defined  through  specification  of 

•  The  number  of  components 

•  The  relative  volume  of  each  component 

•  The  yield  strength  for  each  component 
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Fig.  3-3  Stress  Strain  Curve  Showing  Bauschinger  Effect 


(2)  The  strains  are  estimated  for  all  points  in  the  shell  over  the  shell 
coordinates  and  through  the  thickness.  This  generally  is  done  through 
extrapolation  from  previous  solutions. 

(3)  A  subroutine  is  called  within  which,  for  each  of  the  material  components, 
the  stress  corresponding  to  the  assumed  strains  is  determined.  The  total 
stress  for  the  composite  then  is  found. 

(4)  Once  total  strains  and  stresses  are  known,  the  plastic  part  of  the  strain 
increment  can  be  determined  and  added  as  a  pseudo-load  in  an  elastic 
analysis. 

(5)  New  strains  are  computed  and  used  as  estimates.  The  procedure  is  con¬ 
tinued  until  the  computed  strains  agree  to  within  a  given  margin  with  the 
estimated  strains. 


The  following  operations  are  performed  in  the  above-referenced  subroutine: 


(1)  Information  about  material  properties  is  transferred  into  the  routine  to¬ 
gether  with  the  estimated  strain  increments  ( Af  j ,  A^  »  and  Ay)  and 
stresses  at  the  end  of  the  previous  load  step  (o.  ,  j  ,  r  ). 

1  u 

(2)  New  stresses  are  computed  under  the  assumption  that  the  load  step  is 
elastic. 

+  7T7  (Aei +  1'A£2) 

"2  =S2+777<A£2  +  ,,A£1> 
t  -  r  +  AyE/12  (1  +  v)l 

(3)  Set  e£  ■  -  af,  *  kV 

where  k  is  the  ellipse  ratio  for  the  assumed  yield  surface  (usually  /3 ). 

2  2 

(4)  If  <rT  is  less  than  c r^,  the  load  step  is  elastic  in  this  component  (loading 

or  unloading).  If  this  is  the  case  for  all  components,  the  calculations  for 

the  load  step  are  concluded.  There  are  no  pseudo  loads  caused  by  plastic 

strain  increments. 

2  2 

(5)  If  oT  is  larger  than  try  for  some  component,  the  step  is  at  least  partly 
inelastic  for  this  component.  As  we  have  assumed  ideal  plasticity,  the 
stresses  can  be  computed  from  the  conditions  that 

+  <r\  -  *f2  +  k2r2  =  <7y  (26) 


(24) 


(25) 
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where 


°i  =  ~ai +  -777  l*i  -  +  *  (^2  -  H>>] 

ff2  e  ^2  +  777  [^2  -  ^2  +  "(^1  '  *?>]  (27) 

T  =  f+  2TTTTT  lAy  ‘  AyPl 

and  that  the  plastic  flow  is  perpendicular  to  the  yield  surface 

Ac1}  2cr  -  a  AfP  2a  -  a 

*  _1_2  .  .  _J_2  (28) 

AfP  2a2  -  a2  AyP  2lr  t 

After  the  stresses  have  been  determined  in  the  components,  the  average 
stress  in  the  composite  is  found  readily.  As  the  elastic  constants  are  the 
same  for  all  components,  the  plastic  part  of  the  strain  increment  (i.e. , 
the  pseudo  loads),  can  easily  be  obtained. 

3.  5  FINITE  DIFFERENCE  APPROXIMATIONS 

The  STAGS  program  includes  a  capability  of  using  grids  with  variable  spacing.  The 
discretization  of  functions  and  their  derivatives  in  such  nets  has  been  handled  in  the 
following  manner. 

Suppose  a  shell  panel  has  been  covered  with  a  system  of  mesh  lines  whose  coordinates 
are  given  by 


x.  ,  i  =  1 ,  m 


and 


(29) 


,  J  =  1 ,  n 
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where  x  and  6  are  the  axial  and  circumferential  coordinates,  respectively.  Corre¬ 
sponding  to  each  pair  of  values  (i ,  j )  a  rectangular  region  R.  .  is  defined  with  sides 

*  i  J 

of  length 


1/2 


'i  +  1 


-  x. 


i-1 


■  w  iVi 


(30) 


The  regions  R.  .  (and  lengths  a.  .  ,  b.  . )  are  modified  at  boundaries  of  a  shell 
*  *  J  1 1  J  1 1  ] 

by  considering  only  those  portions  which  would  be  within  the  panel.  A  double  integral 
of  a  function  f  over  the  region  R  of  the  panel  may  then  be  approximated  by 


// 


m  n 

-  i  i 

i  =  1  j  =1 


(31) 


The  discretization  is  completed  when  the  integrand  functions  f.  .  are  evaluated  at 

*  *  J 

the  centroids  of  the  regions  R^  ^  in  terms  of  the'neighboring  displacement  components. 

First,  note  that  the  tangential  displacements  u  and  v  have  been  located  at  corners 
of  the  regions  R^  j.  Furthermore,  the  energy  expressions  for  a  general  shell  in¬ 
clude  derivatives  of  u  and  v  only  up  to  the  first  order.  Hence,  even  with  arbitrary 
rectangular  spacing,  only  central  difference  formulas  for  the  u  and  v  displacements 
are  required.  In  contrast,  the  normal  displacement  w  has  been  located  at  the  mesh 
node  points  (x^  ,  6^ ),  and  more  general  finite  difference  formulas  must  be  developed. 

The  coordinates  of  the  centroid  of  a  region  R.  .  are  given  by 

1 1  J 

x.  =  1/4  (Xj_j  +  2x.  +  xI  +  1) 

3j  =  1/4  (flj.j  +  20.  +  0j+1)  (32) 
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Variable  spacing  is  considered  first  with  respect  to  a  single  coordinate  x  only.  With 
the  help  of  a  Taylor's  expansion  (or  equivalently  by  the  use  of  interpolation  formulas), 
the  difference  formulas  for  w  ,  w ,  and  w  >xx  at  x^  may  be  established  as 

(w)i  s  wlx.  =  wl-l/16  ’  [(h  "  k)  *  (3k  +  h)/(hZ  +  hk)1 

+  w^ie  •  ( (h  +  3k)  •  (3h  +  k )/h  +  k))  (33) 

+  wi+1/16  •  [<k  -  h)  •  (3h  +  k )/( hk  +  k2)] 

<w*x>i  a  w'xlxi  "  '  Wi-l/(2h) 

+  (l/(2h)  -  l/( 2k ) ]  (34) 

+  wi+  j/(2k) 

(w’xx>i  s  w'xxlXi  ■  wi-l  •  2/Ih  *  <h  +  k>l 

-  Wj  •  2/(h  •  k)  (35) 

+  w1+1  •  2/(k  •  (h  +  k)J 


where 


h 


xi  ‘  xi-l 


k 


xi+ 1  "  xi 


(36) 


The  corresponding  formulas  for  the  0  coordinate  are  obtained  by  appropriate  substi¬ 
tutions  and  are  denoted  with  superscripts 


(w, 


(w. 


e  ^  ~  w  •  e  1 3 j 
ee^  H  w ’  ee I  3 


(37) 
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im-' 


The  required  two-dimensional  difference  formulas  now  are  obtained  by  combining  the 
formulas  for  both  coordinate  directions 


*i.j  'wl(x..5.)  -  «w)i>)  ’  ««A 


W  ,  XX.  ,=W,  I  =((w,  ))' 

‘•i  “V 


J 


(38) 


W'091,J  ”  w  ■  os  I 


W,x8.  .  5  w,  .1  =  ((w,  )  )■ 

l’J  X0  (X.,  S{)  Xi ,9 


) 


In  general ,  these  equations  involve  the  9  point  "star"  of  neighboring 
values.  However,  it  is  easily  seen  that  all  of  the  formulas  reduce  to 
the  standard  central  difference  formulas  when  uniform  rectangular 
spacing  is  used.  All  of  the  difference  formulas  are  exact  when  the 
displacement  function  w  behaves  quadratically . 

3.6  SHELL  GEOMETRY 


' - V 

9 

Geometry  routines  are  included  in  the  program  for  the  following  common  shell  types: 


•  Plate 

•  Cylinder 

•  Cone/Annular  Plate 

•  Sphere 

•  Torus 


9  Ellipsoid 

•  Parabolid 

•  Hyperboloid 

•  Elliptic  Cylinder 

•  Elliptic  Cone 


These  are  denoted  by  the  parameter  NSHELL  =  1  through  10  as  discussed  in  Section  u. 
In  addition,  two  dummy  routines  (NSHELL  =11  and  12)  are  included  that  are  used  for 
shell  geometries  not  listed  above.  These  enable  the  user  to  provide  his  own  geometry. 
The  following  are  instructions  for  the  generation  of  such  routines. 
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A  set  of  coordinates,  the  surface  coordinates,  are  defined  so  that  a  pair  of  values  (£,tj) 
uniquely  define  the  position  of  a  point  on  the  shell  surface.  and  the  outward  normal 
n  should  conform  to  the  right-hand  rule.  For  example,  for  a  cylinder  it  is  practical  to 
choose  the  axial  and  the  angular  coordinates.  In  addition,  a  set  of  orthogonal 
Cartesian  coordinates  x,  y,  x  is  defined,  and  the  coordinate  values  in  this  system 
for  a  point  on  the  shall  surface  are  expressed  in  terms  of  surface  coordinates. 
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The  derivatives  of  these  coefficients  with  respect  to  the  surface  coordinates,  4,17  , 
are  also  required.  These  are  written 


AX  =  dA/di  ,  BX  =  9B/3|  . FX  =  9F/3£ 

AY  =  aA/arj,  BY  =  9B/87),...,  FY  =  dF/dT , 


(44) 


If  the  surface  coordinates  |,tj  are  nonorthogonal,  NSHELL  is  set  equal  to  12  and  the 
user  provides  a  subroutine  UNORTH  (see  Section  5.5).  if  the  coordinates  are  ortho¬ 
gonal,  the  procedure  is  simplified  since  C  =  E  =  0.  In  that  case  NSHELL  should  be 
set  equal  to  11  and  the  user  provides  a  subroutine  ORTH  (see  Section  5.4).  Examples 
of  the  subroutines  ORTH  and  UNORTH  are  presented  below. 


Consider  the  development  of  the  subroutine  ORTH  for  a  paraboloidal  shell  with  meridian 
defined  by 


x  =  ar2  (45) 

A  cylindrical  system  with  the  axial  coordinate  (  ahd  the  angular  coordinate  r\  can 
be  chosen  as  surface  coordinates.  The  Cartesian  coordinates  are  expressed  in  terms 
of  the  surface  coordinates: 


x  =  k  .  V )  =  ( 

y  =  g(*.’?)  =  U/a)1/2  sin  ij  (46) 

z  =  h( ?  ,  tj )  =  (4/a )1//2  cos  tj 

Hence 


1 

Si 

=  0 

sin  ij/(  2 

g.  u 

=  -  sin  r\/ ( 4| 

cos  tj/( 2  \/7a") 

h,  ^ 

=  -  cos  »)/(44  / 4a) 
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0 

*-n 

-3 

-3 

« 

O 

(4/a)1/2  cos 

s.  „„  -  -  (?/a)1/28i«n 

-  (4/a)V2 

1/2 

h,  »  -  '4/a)  cos  7) 

7)7) 

(47) 

The  components  of  the  normal  are 

nx  =-(sin  rj/(2  /a4 )  •  sin  r?  4J/1 

+  cos  tj/(  2  4 a| )  •  cos  i )  4 f/a)^g 

n2  =  (sin  r,  •  /T7i)  ^ 
n3  =  (cos  77  •  /!7a)3S 

Substitution  of  equations  (47  and  48)  into  (43a-f,44)  yields 


0.5  _1 
a  .B 


(48) 


D 

F 

DX 

DY 


*  t  .  leilLJl  +  Icos 2_jl 

A  l  4  4  |a 


1/2 


*/l  +  l/(4a4  ) 


B  ~  [  a  cos^  ^  +  a  = 

c  =  0  /  ,  \ 

AX  =  |  [  1  +  l/(4a4  )] 1/2  ^  (-  p  j  =  -  1/144  Va4  (1  +  4a4  )) 


(49) 


BX  =  1/(2  /af) 

AY  =  BY  =  0  =  CX  c  CY 


1 

AB 


{W-  r3/2-l-2  ,-i-i)--V(4aABC 


t  \ 


^  (-  •** » ■  4  -  °os>  -I)*-  {/<*AB' 


(50) 


1  I/IX-B  *  BX-  AH  4  AB  FX  =  .  >  A^-.JAX-B^  BX-AK. 
'  “  - ^V2  “  <“>' 


4a 

=  0 


FY  =  0 
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Table  4  (Section  5)  shows  the  corresponding  subroutine  ORTH  for  this  example. 

As  an  example  of  a  shell  with  nonorthogonal  coordinates,  the  development  of  the  Sub¬ 
routine  L’NOUTH  for  an  elliptic  ll  cone  is  given  here. 

The  tollowing  sketch  provides  the  geometry  of  an  elliptical  cone.  The  shell  surface 
and  its  boundaries  are  defined  by  the  distances  a,  B,  and  c.  The  geometrical  con¬ 
stants  occurring  the  kinematic  relations  can  most  conveniently  be  determined  if 
an  elliptical  coordinate  system  is  used  in  the  definition  of  the  shell  midsurface.  Hence 
with 


h*-—' 


a  =  a/c  semi -major  axis 
b  =  B/c  semi-minor  axis 


the  coordinates  in  a  Cartesian  system  x ,  y ,  z  can  be  written 


'51) 


x  =  a|  cos  Tj 

y  “  b|  sin  t)  (52) 

z  =  f 

Tne  surface  coordinate  £  represents  the  distance  from  the  apex  of  the  cone  and  tj 
can  be  expressed  in  terms  of  the  angular  coordinate  <p  through 
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Tj  ■  arctan  tan  <f>J  (53) 

The  coordinate  lines  intersect  one  another  in  an  angle 

6  «  arc  cos  { C/(  AB )  ]  (54) 

The  coefficients  in  the  first  fundamental  form  becomes 

A  =  ( 1  +  a2  cos2  i)  +  b2  sin2  ij),/2 

C  (a2  -  b2)sin  rjcost)  (55) 

B  ■  4  ( a2  sin2  tj  +  b2  cos2  tj  )^2 

The  derivatives  of  these  coefficients  are 

AX  «  0 
CX  =  c/{ 

BX  «  B/{ 

AY  =  -  ( a2  -  b2 )  sin  tj  cos  rj/A 
CY  «  -  £(a2  -  b2)  (cos2  tj  ~  sin2  t?) 

BY  *  £2  sin  rj  cos  ?j(a2  -  b2  )/B  (56) 

The  coefficients  of  the  second  fundamental  form  and  their  derivatives  are 

D  -  E  =  0 

DX  ®  EX  =  0 
DY  -  EY  =  0 
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F 


+  ab£2/H 


FX  =  F/{ 

FY  =  -ab£  (£/H  )2  sin  tj  cos  tj  ( a2  -  b2 )  (57) 

where 

H  =  Ha2  sin2  tj  +  b2  cos2  r>  +  a2b2  )1/2  (58) 

Table  5  (Section  5)  shows  the  subroutine  UNORTH  for  an  elliptical  cone.  Note:  in  this 
example  the  location  of  columns  in  the  output  will  be  in  r/  coordinate  rather  than  in 
angular  coordinate  <p  • 

3.7  LOADING 

The  program  will  accept  input  defining  two  independent  load  systems.  During  one  run, 
the  loads  in  the  two  systems  may  both  be  proportionally  increased,  or  either  one  can 
be  restrained  by  use  of  the  "Maximum  Load"  input.  The  collapse  analysis  for  a  shell 
may  consist  of  several  runs,  the  restarting  capability  is  an  essential  feature  of  the 
computer  program.  This  restarting  capability  also  allows  the  user  to  change  the  rela¬ 
tion  between  the  increments  of  the  two  load  systems  at  any  given  level  of  loading.  This 
makes  it  (for  example)  possible  to  apply  a  fixed  external  pressure  and,  in  subsequent 
runs,  to  keep  the  pressure  constant  while  an  axial  shortening  is  applied  until  collapse 
occurs.  The  loading  can  be  either  a  prescribed  load  or  a  prescribed  displacement  at 
any  of  the  grid  points.  In  the  bifurcation  analysis,  the  critical  loading  is  defined  as 
the  initial  load  for  system  B  +  eigenvalue  times  the  initial  load  for  system  A.  Con¬ 
sequently  if  we  want  to  find  the  critical  axial  load  of  a  shell  with  a  fixed  internal 
pressure  we  can  represent  the  internal  pressure  by  load  system  B  and  the  axi**i  com¬ 
pression  by  load  system  A. 

The  load  systems  can  contain  uniform  surface  tractions  or  line  loads  in  any  direction. 
The  line  loads  can  be  applied  in  the  form  of  a  specified  load  or  a  specified  displacement 
along  any  of  the  boundary  lines.  Point  loads  or  w  displacements  can  be  defined  at  any 
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grid  point.  In  addition,  the  user  has  the  option  to  add  a  subroutine  which  specifies  non- 
uniform  pressure  or  line  loads  as  functions  of  the  shell  coordinates.  Such  nonuniform 
loads  are  internally  converted  into  point  forces  at  the  grid  points. 

3.8  CONSTITUTIVE  RELATIONS 

The  constitutive  equations  are  here  used  in  the  following  form. 


'N1 ' 

N2 

C2 

N12 

*12 

< 

Mi 

>  -  [€„)  « 

► 

*1 

M2 

*2 

*12. 

.2*12. 

The  stiffness  coefficients,  ,  can  be  computed  from  the  elastic  and  geometric  prop¬ 
erties  of  the  shell  wall.  Special  subroutines  in  which  these  properties  are  read  and 
the  stiffness  coefficients  computed  are  provided  for 

1.  Monocoque  orthotropic  shells 

2.  Shells  with  skew  stiffeners  (waffle  pattern) 

3.  Fiber  reinforce'* 

4.  Layerf'?  *'  _  ^nnotropic  layers) 

"  oaells  with  corrugated  skin 

6.  Shells  with  one  corrugated  and  one  smooth  skin 

7.  Shells  in  which  the  elastic  properties  vary  gradually  through  the  thickness 

In  addition,  a  subroutine  is  provided  in  which  the  stiffness  coefficients  are  modified  to 
Include  the  effects  of  "smeared"  stiffeners. 

For  generality,  an  option  is  provided  in  which  the  user  writes  his  own  wall  property 
routine  (WALL).  Instructions  are  given  here  for  the  derivation  of  such  subroutines; 
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subsequent  text  provides  an  example.  This  option  provides  the  only  method  for  intro¬ 
ducing  variable  wall  properties  (except  variable  elastic  properties  in  monocoque  shells; 
see  input  pr  a  rat  ion  section  under  card  M-l  for  IWALL  -  8  or  9).  The  derivation  of 
the  user-written  wall  property  routine  (WALL),  has  been  facilitated  because  standard 
subroutines  for  shell  types  listed  above  can  be  called  from  WALL. 

For  example,  in  his  routine  the  user  can  define  such  properties  as  modulii,  shell 
thicknesses,  and  stiffener  data  as  functions  of  the  surface  coordinate,  and  call  the 
appropriate  shell  wall  subroutine  which  will  return  the  corresponding  stiffness  coef¬ 
ficients.  The  ’.iser-written  subroutine  will  automatically  be  called  whenever  stiffness 
coefficients  are  needed  if  the  parameter  IWALL  (M-l  card)  is  set  equal  to  one.  More 
than  one  of  the  special  routines  may  be  called  from  the  user  written  routine  if  the  shell 
wall  Is  of  different  type  in  different  areas. 

For  a  case  in  which  the  shell  wall  is  of  a  type  not  included  among  the  standard  types, 
stiffness  coefficients  (C^  matrix)  may  be  read  directly  in  WALL.  The  coefficients  may 
be  derived  as  described  in  the  following.  The  strain  energy,  U,  is  expressed  in  terms 
of  strains  and  changes  of  curvature.  If  a  subscript  following  a  comma  indicates  differ¬ 
entiation  with  respect  to  one  of  the  strains  and  changes  of  curvature  such  that  for 
i  =  1,  2,  3,4, 5, 6  derivatives  are  taken  with  respect  to  e  \ »  ’  ^12  ’  *1  ’  K2  ’  ^K12  ' 

respectively,  then 


(59) 


Stiffness  coefficients  for  the  standard  wall  types  are  given  in  Ref.  14  which  also  gives 
as  an  example  the  detailed  derivation  for  one  case,  the  fiber  reinforced  shell.  As  a 
somewhat  simpler  example  we  will  demonstrate  here,  how  the  coefficients  are  obtained 
for  a  shell  with  rectangular  stiffeners  in  one  direction.  It  is  assumed  that  the  effects 
of  torsional  stiffness  and  resistance  to  rotation  of  the  stiffener  can  be  neglected. 


If 


e  ,  e  ,  and  y 
x  ’  y  ’  r\y 
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denote  strains  at  the  reference  surface  (here  the  midsurface  of  the  skin)  and 


are  the  changes  of  curvature, 
can  be  obtained  from 


k  ,  k  ,  and  k 
x  y  xy 


then  the  strains  at  any  point  off  the  reference  surface 
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X  X 

€+  k  Z  in  skin 

y  y 
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irrelevant  in  stringer 


The  geometric  properties  of  the  shell  wall  are  shown  in  Fig.  3-4. 


Fig.  3-4 

The  strain  energy  per  unit  area  of  the  shell  wall  is 

U  =  ; A  /  (a  e  +  a  €  .  +  t  V  ^dv  (60) 

2d  J  \  x  x  y  y  xy  rxy/  '  ' 
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The  contribution  from  the  skin  is 
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For  the  stringer  we  have 
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or  with 
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=  bh  =  stringer  area 
„  bh2/12  =  stringer  moment  of  inertta 

-  (h  +  t)/2  =  stringer  eccentricity 


=  —  [ 
USTR  2d  [' 


Ai2  +  2Ae€_xKx  + 


2  2l 

(I  +  Ae>l 


(63) 


Then  by  use  of  Eq.  (59)  we  find 
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If  a  user  written  shell  wall  property  routine  is  used, certain  rules  must  be  followed. 
Thus  the  list  of  variables  under  the  subroutine  name  (WALL)  is  (X,Y,CCC)  which 
represent  the  shell  surface  coordinates  and  stiffness  coefficients.  The  X  and  Y 
dimensions  must  correspond  to  those  shown  in  Section  5,  page  5-1. 


If  any  of  the  other  routines  are  called  from  WALL  the  list  in  the  call  statement  must  be 
(N) .  Any  number  of  the  routines  listed  In  Table  7  may  be  called  from  WALL.  The 
table  also  shows  the  required  common  statements  corresponding  to  each  of  the  sub¬ 
routines.  As  there  are  no  cards  to  be  read  in  the  routines  called  from  WALL  it  is 
necessary  to  provide  in  WALL  the  size  of  the  quantities  in  these  common  statements. 
They  are  defined  on  the  cards  M-2B  through  N-4B  in  the  input  description.  If 
smeared  stiffeners  data  are  provided  in  WALL  the  user  must  CALL  STIFF  in 
subroutine  WALL  to  generate  the  stiffeners  coefficients. 


As  an  example  we  will  consider  here  a  shell  which  is  composed  of  a  spherical  segment 
and  a  cylindrical  part  as  shown  in  Fig.  3-5. 
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AXIS  OF 


The  user  written  geometry  which  is  needed  for  this  case  we  will  assume  uses  as  X 
coordinate  the  arc-length  from  the  upper  edge.  Hence  the  local  radius 

r  =  10  •  cos  (J  -  i) 

if  there  are  40  stringers  around  the  circumference  we  have 

,  2  r  •  r  _  7r  •  r 

40  20 

We  assume  that  it  is  specified  that  the  stringer  height  varies  linearly  from  zero  at  the 
shell  edge  to  0.  6  inch  at  the  interface  between  the  spherical  and  cylindrical  shells. 


3-30 


h 


X‘0.6  _  2.4  •  X 
(TT/4)  “  7T 

7 

The  material  for  both  parts  of  the  shell  is  isotropic  with  E  =  10  and  v  -  0.3  . 
The  routine  WALL  corresponding  to  this  shell  is  listed  as  an  example  in  Table  5-7. 
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Section  4 
STRATEGY 


The  preceding  paragraphs  have  described  the  scope  and  limitations  of  the  STAGS 
computer  program.  The  following  text  discusses  the  choice  of  certain  control  param¬ 
eters;  and  an  effort  is  made  to  convey  to  the  user  some  of  the  experience  that  has 
been  gained  through  extensive  program  use. 

There  are  essentially  two  areas  within  which  the  user  is  required  to  use  his  own 
judgment.  The  first  of  these  is  the  modeling  of  a  structure  so  that  it  becomes  ame¬ 
nable  for  analysis.  This  will  include  the  choice  of  a  able  grid.  The  second  area 
is  the  choice  of  certain  parameters,  such  as  the  initial  load  step  and  the  convergence 
criterion,  which  will  govern  the  flow  of  computations  i/ie  nonlinear  analysis. 

The  modeling  of  the  structure  is  generally  consider  be  outside  of  the  scope  of  this 

manual.  It  is  assumed  that  the  user  has  decided  abo’  the  shell  geometry  and  material 
properties,  about  boundary  conditions  and  loading,  and  about  what  to  include  in  terms 
of  stiffeners  and  cutouts.  The  next  step  would  be  to  determine  which  type  of  analysis 
to  apply.  If  the  interest  is  in  the  stress  distribution  for  loads  which  are  known  to  be 
small  in  comparison  to  the  collapse  load  of  the  shell  a  purely  linear  analysis  can  be 
used.  If  it  can  be  safely  assumed  that  changes  in  geometry  or  stress  distribution  are 
negligible  at  loads  only  slightly  below  collapse,  the  branch  for  bifurcation  buckling  can 
be  used  for  establishment  of  the  stability  limit. 

In  case  there  is  some  uncertainty  about  this,  the  following  procedure  might  help  the 
user  to  avoid  an  expensive  complete  nonlinear  analysis.  The  design  load  is  used  as 
both  initial  and  maximum  load.  The  linear  analysis  then  will  be  obtained,  and  an 
attempt  will  be  made  at  solution  of  the  nonlinear  equations  at  this  load  step.  If  the 
linear  analysis  represents  a  good  approximation,  convergence  will  be  obtained  within 
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a  few  Iterations  and,  as  no  additional  factoring  is  required,  this  will  add  very  little  to 
the  computer  time.  If  convergence  is  not  obtained,  this  will  serve  as  a  warning  that  a 
nonlinear  analysis  may  be  needed.  It  is  recommended,  therefore,  that  the  nonlinear 
option  be  used  unless  similar  cases  previously  analyzed  clearly  indicate  that  a  linear 
analysis  is  satisfactory.  Also,  note  that  a  slight  increase  in  load  beyond  the  design 
load  may  lead  to  a  considerable  increase  in  the  influence  of  nonlinear  terms.  Such  be¬ 
havior  would  be  revealed  by  a  bifurcation  analysis.  If  results  from  the  nonlinear  analy¬ 
sis  at  the  design  lead  differ  little  from  results  from  a  linear  analysis  and  the  bifurcation 
load  is  well  above  the  design  load,  shell  collapse  is  not  a  problem.  In  choosing  the  type 
of  analysis,  it  may  be  useful  for  the  program  user  to  know  that  a  nonlinear  analysis,  if 
it  converges  at  the  first  step  only  requires  about  25  percent  more  computer  time  than 
a  linear  analysis,  and  that  a  bifurcation  buckling  analysis  could  take  from  two  to  three 
times  as  much  computer  time  as  a  linear  analysis,  when  multiple  eigenvalues  exist  in 
the  neighborhood  of  the  lowest  eigenvalue . 

If  the  shell  is  subjected  to  such  load  as  to  make  instability  a  possibility,  and  if  the 
shell  dimensions  are  not  determined  by  other  -  nonstructural  -  requirements,  efficient 
design  can  be  achieved  only  if  a  collapse  analysis  is  performed.  In  this  case,  the  bi¬ 
furcation  buckling  analysis  should  be  used  only  if  the  analysis  of  similar  cases  clearly 
indicates  its  adequacy. 

The  next  step  is  to  determine  the  size  of  the  grid.  For  economy,  it  is  necessary  to 
take  advantage,  whenever  possible,  of  the  capability  of  using  nets  with  variable  spacing. 

If  similar  cases  have  not  been  run  before,  it  is  essential  that  a  convergence  study  be 
made.  Often  this  may  be  done  by  use  of  a  linear  analysis  or  a  bifurcation  buckling 
analysis.  Note,  however,  that  sometimes  a  good  estimate  of  the  collapse  load  can  be 
obtained  with  e  more  coarse  grid  than  is  required  for  an  accurate  estimate  of  pre¬ 
buckling  stress  distribution.  Also,  sometimes,  the  coarseness  of  a  net  which  is  satis¬ 
factory  for  linear  analysis,  will  result  in  spurious  buckling  modes. 

The  user  can  reasonably  well  determine  the  computer  time  for  a  linear  stress  analysis 
or  bifurcation  buckling  analysis  from  Fig.  4-1.  For  a  nonlinear  analysis,  it  is  difficult 
to  predict  how  many  refactorings  are  needed  and  how  big  the  local  steps  may  be. 
Depending  upon  how  drastic  the  changes  in  geometry  or  stress  distribution  are  in  the 
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RUN  TIME  (SEC) 


Fig.  4-1  CDC  6600  Computer  Time  for  Bifurcation  and  Linear  Stress  Analysis 
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critical  range,  the  computer  time  for  a  collapse  analysis  varies  from,  for  example, 

10  to  100  times  the  time  for  a  linear  analysis. 

Once  the  structural  model  (including  the  grid)  has  been  chosen,  the  linear  analysis  is 
straightforward;  for  a  nonlinear  analysis  additional  decisions  must  be  made.  The  user 
must  choose  the  initial  load  step  and  other  control  parameters  and  also  may  elect  to 
override  built-in  values  for  the  iteration  convergence  criterion  and  the  over-  or 
under-relaxation  factor.  An  effort  has  been  made  to  make  automatic  as  much  as  pos¬ 
sible  of  the  computation  strategy.  Thus,  for  example,  if  convergence  is  slow,  an 
over-  or  under-relaxation  factor  is  automatically  applied  depending  on  whether  the 
convergence  is  uniform  or  oscillating.  The  user  generally  leaves  blank  fields  for 
convergence  parameter  and  relaxation  factor.  This  means  that  the  convergence 

_4 

parameter  is  set  to  10  ,  and  the  relaxation  parameter  is  chosen  as  indicated  above. 
The  user  may  loosen  this  convergence  criterion  if  .t  is  indicated  from  the  results  that 
round-off  errors  make  convergence  difficult  or  ever,  impossible.  The  criterion  also 
may  be  tightened  as  discussed  below.  The  relaxation  factor  should  be  determined  by 
the  user  only  if  the  user  has  considerable  experience  in  the  use  of  the  program.  If  the 
relaxation  factor  is  given  any  value  on  the  input  carB  (including  1.  0),  this  factor  will  be 
used  throughout  the  analysis. 

The  linear  solution  is  always  computed  and  used  as  a  first  estimate  for  the  unknowns. 
Later,  estimates  are  obtained  through  extrapolation  (quadratic  when  two  or  more  pre¬ 
vious  solutions  are  available).  Therefore,  the  initial  load  must  be  chosen  so  that  the 
solution  for  this  load  level  is  not  too  far  from  the  linear  solution.  A  bifurcation  analy¬ 
sis  to  guide  in  the  choice  of  initial  stepslze,  although  expensive,  is  possible  if  there  is 
no  other  basis  for  a  reasonable  guess.  If  convergence  is  obtained  easily  within  one  or 
two  Iterations,  the  load-step  will  be  automatically  increased  by  a  factor  of  1. 5.  This 
increase  will  be  repeated  any  number  of  times,  and  thus  the  penalty  for  a  poor  initial 
load  step  guess  will  be  alleviated.  This  feature  in  the  program  also  is  useful  for  many 
cases  in  which  a  proper  guess  is  made  for  starting  values.  For  example,  for  a  shell 
of  revolution  with  nonuniform  load,  convergence  is  usually  slow  at  the  first  couple  of 
load  steps;  after  previous  solutions  are  available  for  extrapolation,  the  convergence  is 
faster  and  the  load  step  can  be  increased. 
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The  choice  of  initial  load  and  of  initial  size  of  the  load  step  for  maximum  efficiency  is 
case  dependent.  A  user  who  is  experienced  in  nonlinear  analysis  can  choose  these 
parameters  to  suit  his  particular  case.  To  a  user  who  has  no  experience  on  which  to 
base  such  a  decision,  one  might  recommend  that  initial  load  and  initial  step  be  chosen 
at  about  one  tenth  of  the  anticipated  collapse  load.  If  the  structure  is  expected  to  behave 
almost  linearly  until  the  critical  load  is  closely  approached,  the  initial  load  may  be 
large  and  the  initial  step  much  smaller. 

If  convergence  is  not  obtained  in  the  first  load  step,  it  is  probably  because  the  initial 
loads  are  too  big,  and  the  nonlinear  solution  differs  too  much  from  the  linear  solution. 
However,  there  are  other  possibilities  such  as  simple  mistakes  in  input  quantities. 

The  user  also  is  advised  to  check  whether  the  boundary  conditions  allow  rigid  body 
motions  -  translations,  rotations,  or  combinations  of  these.  It  also  is  poP'dble  that 
the  boundary  conditions  allow  deformation  in  a  mode  in  which  the  shell  is  very  weak, 
e.g.  an  inextensional  deformation  mode.  For  such  cases,  it  is  possible  that  the  sys¬ 
tem  is  so  ill  conditioned  that  the  problem  cannot  be  solved  with  single  precision  accu¬ 
racy.  If  this  deformation  mode  is  not  essential  to  the  solution  of  the  problem,  it  can 
be  restrained  and  subsequently  a  successful  analysis  can  be  executed.  The  >'ases  in 
which  such  difficulties  occur  are  quite  rare. 

If  the  convergence  was  not  obtained  because  the  initial  load  step  was  too  large,  then  this 
step  must  be  reduced.  If  the  Iteration  diverges,  the  load  may  be  cut  by  as  much  as  a 
factor  of  10,  whereas  a  lesser  reduction  may  be  successful  if  the  iteration  converges 
too  slowly.  After  initial  convergence  has  been  obtained,  the  program  continues  with 
the  input  load  step  until  a  load  level  is  reached  at  which  the  program  fails  to  converge 
within  prescribed  limits.  The  program  then  either  cuts  the  load  increment  in  half  or 
makes  a  standard  Newton  iteration  step  by  recomputing  L'(Xm)  where  Xm  is  the  dis¬ 
placement  vector  solution  obtained  for  the  previous  load  step.  The  choice  between  these 
two  alternatives  is  based  on  strategy  parameters  input  by  the  user: 

ICUT  -  Total  number  of  times  load  increment  may  be  cut  in  half 

INEWT  -  Number  of  times  the  matrix  L'(Xra  )  may  be  recomputed 

ISTRAT  -  Number  of  times  load  increment  is  cut  between  recomputation  of 

L,(Xm>- 
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Thus,  if  ISTRAT  =  1,  the  program  alternates  between  cutting  the  load  increment  and 
recomputing  the  matrix  when  convergence  difficulties  are  encountered.  When  the  pro¬ 
gram  has  already  cut  the  step  size  ICUT  times,  additional  "Newton"  steps  (recomputa¬ 
tion  and  factorization  of  L'(  Xm  ) )  are  permitted  until  the  total  exceeds  INEWT.  Of 
course  "he  reverse  situation  also  may  occur  in  which  the  program  is  only  permitted 
to  cut  the  step  size  but  not  to  take  additional  "Newton"  steps.  If  neither  strategy  is 
allowed,  the  program  stops  as  soon  as  convergence  difficulties  arise. 

In  view  of  the  highly  unpredictable  nature  of  nonlinear  behavior,  it  is  very  difficult  to 
prescribe  the  oest  values  for  the  strategy  parameters  in  advance.  As  an  initial  choice 
for  an  unfamiliar  type  of  shell,  the  values 

ICUT=2 
INEWT  =2 
ISTRAT-1 

are  suggested.  By  observing  the  convergence  behavior  in  previous  computer  runs,  a 
more  effective  set  of  values  may  be  selected  for  continuation.  Thus,  if  for  example, 
the  last  five  load  steps  converge  in  one  iteration  each,  a  larger  load  increment  may  be 
desirable.  Conversely,  as  a  collapse  load  is  approached,  it  will  be  necessary  to  cut 
the  load  increment  frequently  (e.g. ,  ISTRAT=2  or  ISTRAT=3).  Sometimes  a  shell 
exhibits  substantial  redistribution  of  stresses  at  loads  well  below  the  collapse  load. 

In  such  regions,  the  most  economical  strategy  may  be  to  take  "Newton"  steps  (which 
take  into  account  the  stress  redistribution)  more  frequently  than  cutting  the  step  size. 
Such  **.  strategy  could  be  achieved  by 

ICUT=1 
INEWT  *6 
ISTRAT=1 

In  general,  it  has  been  found  advisable  to  restrict  each  computer  run  to  5  to  10  times 
the  computer  time  for  factoring  so  that  the  convergence  behavior  can  guide  the  selec¬ 
tion  of  strategy. 
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The  STAGS  program  is,  of  cours?.  not  an  appropriate  tool  for  analysis  of  the  buckling 
or  collapse  of  shells  Oi  revolution  under  axi  symmetric  loading.  For  such  problems, 
simpler  tools  are  available,  such  as  the  B0S0R4  program  (Ref.  15).  Axisymmetric 
cases  have  been  analyzed  for  verification  of  the  validity  of  the  STAGS  program. 
Theoretically,  as  the  critical  load  is  reached,  the  round-off  errors  should  trigger  a 
deformation  in  the  buckling  mode.  In  practical  application,  it  is  generally  found  that 
in  the  early  stage  of  buckling  the  amplitude  of  the  buckling  mode  is  so  small  in  com¬ 
parison  to  the  prebuckling  displacement  that  its  growth,  although  in  a  relative  sense 
large,  will  not  violate  the  specified  convergence  criterion.  This  difficulty  is  avoided 
by  the  specification  of  initial  imperfections  in  the  shell  geometry  which  are  small 
enough  not  to  appreciably  affect  the  buckling  load  but  large  enough  to  trigger  the  new 
deformation  pattern.  The  program,  therefore,  has  been  equipped  with  an  option  to 
add  a  subroutine  which  describes  an  initial  lateral  displacement  pattern. 

It  has  been  found  during  use  of  the  program  that  the  imperfections  may  be  useful  as 
triggers  also  in  other  cases  than  those  with  perfect  axial  symmetry.  For  example, 
in  the  analysis  of  an  elliptical  cylinder  with  an  aspect  ratio  of  1. 5,  it  was  found  that 

-5 

without  trigger  it  was  necessary  to  use  a  very  severe  convergence  criterion  (e  =  10  ) 
but  if  a  small  imperfection  is  added,  the  same  collapse  load  may  be  computed  in  about 
half  the  run  time  with  a  less  severe  convergence  criterion  (e  =  10“^).  If  the  elliptical 
cylinder  has  a  significantly  smaller  aspect  ratio,  it  is  likely  that  the  amplitude  of  the 
buckling  pattern  which  is  present  in  the  prebuckling  displacement  is  too  small  to  act 
as  a  trigger.  In  this  case,  the  computation  of  a  collapse  load  must  include  the  use  of 
a  small  imperfection.  The  choice  of  imperfection  mode  may  often  be  aided  by  knowl¬ 
edge  of  the  bifurcation  buckling  mode. 

If  difficulties  like  these  do  occur,  they  will  be  discovered  when  attempted  refactoring 
at  a  load  level  above  the  collapse  load  leads  to  a  coefficient  matrix  for  the  linear  sys¬ 
tem  which  is  not  positive  definite.  In  such  a  case,  the  user  must  either  sharpen  his 
convergence  criterion  or  introduce  an  imperfection. 
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The  run  may  often  be  saved  if  a  new  run  is  restarted  from  an  earlier  solution;  i.  e. , 
ISTART  is  chosen  to  be  1  or  2  rather  than  3,  which  is  the  value  chosen  under  normal 
conditions.  There  also  may  be  other  reasons  to  suspect  that  an  inaccurate  solution  has 
been  accepted  in  which  case  restart  from  file  1  or  file  2  on  the  data  tape  Is  advisable. 

Note  that  through  use  of  additional  analysis  with  various  degree  of  imperfections,  the 
user  of  the  program  can  get  some  notion  of  the  degree  of  imperfection  sensitivity  of  the 
collapse  load. 

The  input  card  with  strategy  parameters  (P-1B  card)  also  includes  a  parameter  ISEC. 
Occasionally,  during  computations,  a  check  is  made  of  whether  the  elapsed  computer 
time  exceeds  ISEC;  in  this  case,  intermediate  results  are  saved  on  data  tape.  ISEC 
should  first  be  chosen  to  be  a  minute  or  so  less  than  the  time  estimate  at  which  the 
operator  aborts  the  run.  Before  each  refactoring,  the  program  also  checks  that  suf¬ 
ficient  time  is  available  to  make  refactoring  meaningful.  To  refactor  at  the  end  of  a  run 
would  be  wasteful  because  a  restarted  run  begins  with  factoring. 

In  the  bifurcation  analysis ,  inverse  power  iteration  is  used  to  obtain  the  critical  value 
closest  (in  absolute  value)  to  the  initial  shift  point.  The  rate  of  convergence  to  the 
critical  mode  may  be  very  slow  unless  a  shift  of  the  eigenvalue  spectrum  is  used. 

If  the  parameter  ISHIFT  is  set  greater  than  zero,  the  program  automatically  performs 
up  to  a  maximum  of  ISHIFT  eigenvalue  shifts  to  expedite  convergence.  It  is  proposed 
that  if  the  user  has  no  special  reason  to  do  otherwise,  ISHIFT  is  set  to  2  and  ITERAT 
(the  maximum  number  of  Iterations  between  shifts)  is  set  to  20. 

Whenever  the  critical  load  is  reasonably  well  known,  it  is  probably  desirable  to  use  a 
value  somewhat  below  the  expected  buckling  load  as  initial  shift.  In  general,  it  should 
be  noted  that  increased  convergence  rates  are  obtained  at  the  cost  of  a  complete  matrix 
factorization  for  each  eigenvalue  shift.  An  initial  shift  may  sometimes  be  necessary 
to  obtain  a  physically  meaningful  critical  value  which  does  not  correspond  to  the  lowest 
mathematical  eigenvalue.  For  example,  when  a  load  system  results  in  tension  some¬ 
where  in  the  structure,  there  will  usually  be  negative  eigenvalues  which  may  not  be  of 
interest.  Also,  it  may  not  be  convenient  to  eliminate  rigid  body  motion  by  means  of 
boundary  conditions  in  which  case  there  will  be  eigenvalues  approximately  or  exactly 
zero.  In  these  cases,  the  physically  meaningful  buckling  load  can  be  obtained  by  an 
initial  shift  which  is  sufficiently  close  to  the  desired  critical  value. 


4-8 


Section  5 

USER- WRITTEN  SUBROUTINES 


To  extend  the  applicability  of  the  STAGS  computer  program,  the  option  of  several  user- 
written  subroutines  was  provided.  These  subroutines  make  It  possible  for  the  user  to 
communicate  to  the  system  functional  relationships  which  would  be  difficult,  If  not  Im¬ 
possible,  to  define  through  regular  data  card  input.  Do  not  read  input  cards  in  any  of 
these  routines. 

Some  of  the  information  defined  in  the  user-written  subroutines  may,  in  effect,  over¬ 
ride  data  read  in  on  regular  input  cards,  but  none  of  the  user-written  routines  suppresses 
the  reading  of  any  of  these  cards. 

Instructions  and  examples  of  user-written  subroutines  are  given  here  except  for  the 
geometry  routines  (ORTH  and  UNORTH)  and  routine  WALL  which  are  discussed  in 
Section  3. 


The  coordinates  X  and  Y  used  in  the  user-written  subroutines  must  correspond  to  the 
coordinates  used  for  the  description  of  the  shell  geometry  (NSHELL)  as  follows: 


NSHELL  GEOMETRY  X 


Y 


1  Cylinder 

2  Cone/ Annular  Plate 

3  Plate 

4  Sphere 

5  Paraboloid 

6  Elliptic  Cylinder 

7  Ellipsoid 

8  Torus 

9  Hyperboloid 

10  Elliptic  Cone 

11  ORTH  1 

12  UNORTH  ) 


Length 

Length 

Length 

Degrees 

Length 

Length 

Degrees 

Degrees 

Degrees 

Length 


Degrees 

Degrees 

Length 

Degrees 

Degrees 

Degrees 

Degrees 

Degrees 

Degrees 

Degrees 


As  specified  by  user  in 
geometry  routine 
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If  trigonometric  terms  are  used  in  a  user-written  subroutine,  the  arguments  must  be 
in  radians.  That  is,  if  X  and/or  Y  are  in  degrees,  they  must  be  converted  to 
radians  in  the  user-written  subroutine  (see  Section  5. 1  and  Table  1  for  example). 

The  portion  of  the  user-written  subroutine  that  must  be  added  to  the  program  file  is 
clearly  marked  with  an  asterisk  (*)  for  particular  examples,  in  Tables  1  to  7.  The 
subroutines  names,  dimensions,  comments,  common  statements,  and  the  return  and 
end  cards  are  part  of  the  permanent  program  file  and  appear  in  it  in  consecutive  order. 

5.1  FUNCTION  WIMP  (K.X.Y) 

This  routine  defines  initial  imperfection,  if  any,  of  the  shell  surface. 

The  computer  program  uses  only  the  first  derivatives  of  the  imperfection  with  respect 
to  the  two  space  variables.  If,  when  the  subroutine  is  entered,  the  parameter  in  the 
list  (K)  is  eqpal  to  2,  the  derivative  with  respect  to  Y  is  requested  (Y  is  zero  on 
boundary  line  4).  Otherwise,  the  subroutine  should  return  only  the  derivative  with 
respect  to  X  (X  is  zero  on  boundary  line  1). 

Table  1  gives  an  example  using  this  routine  for  a  cylinder  with  L  =  10  and  WQ  = 

0. 00001  sin  (Xir/2L)  cos  (6Y) . 

5. 2  SUBROUTINE  USRLD  (X,  Y,  NROW,  NCOL) 

This  subroutine  serves  to  define  a  functional  relationship  between  the  external  loading 
and  the  X,  Y  mesh  coordinates  and  is  called  only  if  LFLG  =  1  on  the  L-l  type  input 
card.  Additional  loads  can  be  defined  by  use  of  the  L-2  cards. 

X(I)  =  1C  coordinate  of  mesh  point  I 

Y  (I)  =  Y  coordinate  of  mesh  point  I 

NROW  =  Number  of  rows 

NCOL  =  Number  of  columns 
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In  the  process  of  coding  the  USRLD  subroutine,  the  user  has  to  define  the  external  load 
P  at  any  or  all  mesh  points  and  then  issue  a  call  statement 

CALL  FORCE  (L,  M,  N,  P,  I) 


where 

L  Row  number  of  mesh  point  where  P  is  acting 
M  =  Column  number  of  mesh  point  where  P  is  acting 
**  =  Direction  of  P 

1  -  Normal  (Z) 

?  -  Tangential  (Y) 

3  -  Tangential  (X) 

P  =  External  load 
I  =  Loadtypa 

-1  -  Displacement 

1  -  Pei  it  force 

2  -  Lin01  .oad  along  rows 

3  -  Line  load  along  columns 

4  -  Pressure  test 

5  -  Line  pressure  load  (uniform  pressure  only) 


The  pos.tive  load  is  applied  in  the  direction  of  positive  displacement. 


Table  2  provides  an  example  using  this  routine  for  internal  pressure  that  varies  along 
the  Y  coordinate  according  to  P  =  10  [cos  (2Y)  +  1). 


5.3  SUBROUTINE  MATER  (X,  Y  ,  IP.TDEG,  EX,  EY,U,  G,  Al,  A2) 

This  subroutine  defines  the  temperature  and  wall  properties  at  every  mesh  point  and 
point  through  the  thickness  in  the  shell  and  is  called  by  the  program  only  if  IWa  LL  =  8 
on  the  M-l  type  input  card. 

X  =  X  coordinate  of  mesh  point 

Y  =  Y  coordinate  of  mesh  point 

IP  =  Number  of  points  across  the  wall  (numbered  inner  to  outer) 

TDEG  =  Wall  temperature 


5-3 


EX  =  Modulus  of  elasticity  in  X  direction 

EY  =  Modulus  of  elasticity  in  Y  direction 

U  =  Poisson's  ratio  (Mxy) 

G  =  Shear  modulus 

A1  =  Coefficient  of  thermal  expansion  in  X  direction 

A2  =  Coefficient  of  thermal  expansion  in  Y  direction 

An  example  using  this  routine  for  material  properties  as  function  of  a  parameter  C  is 
given  in  Table  3,  where 

C  =  1.0  -  0.4  cos  (Y)  for  Y  <  90 

C  =  1.0  for  Y  >90 

5.4  SUBROUTINE  ORTH  (PROP,X,Y) 

This  subroutine  defines  shell  geometries  described  by  a  set  of  orthogonal  surface 
coordinate  lines.  It  must  be  provided  by  the  program  user  when  the  parameter  NS  HELL 
(see  Section  6,  Card  G-2  is  set  equal  to  11.  If  a  3D  plot  of  geometry  (NCHK.GT.0)  is 
requested,  the  user  must  provide  the  orthogonal  cartesian  coordinates  XG,  YG,  and 
ZG  as  shown  in  Table  4. 

An  example  of  a  paraboloid  is  described  in  Section  3.6  and  the  resulting  subroutine 
shown  in  Table  4. 

5.5  SUBROUTINE  UNORTH  (PROP,  X  ,Y) 

This  subroutine  defines  shell  geometries  described  by  a  set  of  nonorthogonal  surface 
coordinate  lines.  It  must  be  provided  by  the  program  user  when  the  parameter  NSHELL 
(see  Section  6,  Card  G-2)  is  set  equal  to  12,  and  can  be  used  only  in  conjunction  of 
IWALL  equal  1  or  2.  (IWALL  =  1  requires  user  written  subroutine  for  shell  wall 
properties. )  If  a  3D  plot  of  geometry  (NCHK.  GT.  0)  is  requested,  the  user  must  pro¬ 
vide  the  orthogonal  cartesian  coordinates,  XG,  YG,  and  ZG  as  shown  in  Table  5. 

Example  of  an  elliptical  cone  is  described  in  Section  3.6  and  the  resulting  subroutine 
is  shown  in  Table  5. 
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5.6  SUBROUTINE  TEMP  (X,  Y,  T,  API,  AP2) 


The  user  bas  the  option  of  entering  the  thermal  loading  in  conjunction  with  load 
Pattern  A  by  means  of  this  subroutine.  Subroutine  TEMP  is  called  for  each  mesh 
point  of  the  shell  and  normally  returns  a  temperature  value  of  zero  unless  the  user 
specifies  otherwise.  X  and  Y  are  the  shell  coordinates  of  the  surface.  API  and 
AP2  are  the  coefficient  of  thermal  expansion  in  X  and  Y  direction,  respectively. 

Table  6  shows  an  example  of  the  use  of  this  routine  for  temperature  variation  given  by 
T  =  50  [1  -  0.4  cos(Y)J  where  Y  is  given  in  degrees. 

5.7  SUBROUTINE  WALL  (X,Y,CCC) 

Subroutine  WALL  makes  it  possible  for  the  user  to  vary  the  stiffness  matrix  or  the 
material  properties  which  influence  the  calculation  of  the  stiffness  matrix  at  various 
mesh-points.  The  common  variables  related  to  the  various  types  of  wall  constructions 
are  available  by  means  of  FORTRAN  COMMON  statements. 

The  subroutine  is  called  at  each  mesh-point  if  IWALL  is  set  to  1  on  the  M-l  type  input 
card.  The  user  may  calculate  the  CCC  stiffness  matrix  (6  x  6)  directly  or  just  set  the 
appropriate  wall  properties  to  the  desired  value  and  call  subroutine  CFB(N)  to  perform 
the  stiffness  matrix  computations  according  to  wall  construction  type  N  (see  Input 
description  for  the  available  wall  construction  options). 

For  a  simple  example,  see  the  description  in  Section  3. 8  and  resulting  subroutine  in 
Table  7. 
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Table  1 

FUNCTION  WIMP 


* 


FUNCTION  W'iNP<)C,X,Y> 

NINP  *  0.0 

P  *  0.00001 

PI2  =  3.14159  /  20.0 

Y1  *  Y  *  3.14159  /  110.0 

IF  <K  ,EQ.  2>  GO  TO  20 

WIMP  «  P  •  PI2  *  COS (PI2  •  X)  *  COS (6.0  •  VI) 
RETURN 

WIMP  ■  -P  *  6.0  •  SIN1PI2  •  X)  •  SINI6.0  •  VI) 

RETURN 

ENO 


Table  2 

SUBROUTINE  USRLD 


* 


l 


SUBROUTINE 
DIMENSION 
00  10  L  * 
DO  10  N  * 
T  *  YIN)  • 
P  «  10.0  * 
CALL  FORCE 


RETURN 

END 


USRLD  (X,  Y,  NROM,  NCOL) 
X(NROW),  Y ( NCOL) 

1,  NROW 

It  NCOL 

3.14159  /  100.0 
( COS <  2 . 3  *  T>  ♦  1.0) 

(L,  M ■  1,  P,  4) 


Table  3 

SUBROUTINE  MATER 

SUBROUTINE  MATER  (X,  Y,  IP,  TQEG,  £X,  EY,  U,  G,  41,  A2I 

OIMENSION  TOEG  < IP) »  EX ( IP) ,EY(IP)  ,UC  I P) , G (IP) , A1 ( IP) ,A2(IP) 

COMMON  FOFST /  TO,Z 

TO  AND  2  MUST  BE  SET  BY  THE  USER 

TO  *  TOTAL  THICKNESS  OF  SMELL  (TD=AT> 

2  »  DISTANCE  FROM  REFERENCE  SURFACE  TO  MIDSURFACE  OF  SMELL  NALL 
f  TO  *  C.l 
2  *  0. 

C  *  1.0 

IF  (Y  .LT.  90.0)  C  ■  1.0  -  0.4  *  COSCY  •  3.14159  /  183.0) 

DO  l  L  *  1,  IP 
TOEG(L)  *  0.0 

<  Exut  *  c  •  icocooo.a 

EY«L)  =  0.4  •  C  •  1000000.0 
U  (L )  *  0.1 

3<L)  *  1.2  •  C  *  1000000.0 

A 1  ( L I  =  0.0 

*2 ( L)  *  0.0 

CONTINUE 

RETURN 

END 
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Table  4 

SUBROUTINE  ORTH 


SUBROUTINE  ORTH  (PROP,  X,  V) 

CIMENSION  PR0P(  3) 

PR0P(1)  *  X  COORDINATE  OF  BOUNDARY  LINE  i 

PROP (2)  *  X  COORDINATE  OF  BOUNOARY  LINE  3 

PROP (3)  =  Y  COORDINATE  OF  BOUNDARY  LINE  4 

PR0P(4I  =  Y  COORDINATE  OF  BOUNDARY  LINE  2 

COMMON  /  FQA  /  NSHELL,  A,  B,  AX,  AY,  BX,  BY ,  XG,  YG,  ZG 
COMMON  /  FQ9  /  C,  CX,  CY,  0,  OX,  DY,  £,  EX,  EY,  F,  FX,  FY 

PARABOLOIO 

»»«»»» »**»»***»»«»**»»«»*9*»» »•»»»»»»»»***»»»» *•»»»»**»»»» 

X  S  s  *  Y  *  Y  EQUATION  OF  PARABOLIC  MERIDIAN 
Ssi • 0 / (4. Q*R)  WHERE  R  IS  THE  DISTANCE  FROM  VERTEX  TO  FOCUS 
PROP ( 5)  *  R  INPUT  ON  CARO  G-2 

PROP(3)  AND  PR0P(4)  ARE  IN  RADIANS  FOR  THIS  EXAMPLE 
S  »  1.0  /  (4.0  *  PR0P(5>> 

X  =  X  ♦  PROP(l) 

Y  =  Y  +  PROP(3) 

XG  -  X 

YG  -  SQRT(X/S)  *  SIN(Y) 

ZG  =  SQRT(X/S)  *  COS(Y) 

T  =  0.25  /  (X  *  S) 

Sfc  A  s  SQRT (1.0  ♦  T) 

B  *  SQRT (x  /  S) 

A1  *  SQRT (S  *  X  *  (1.0  ♦  4.0  *  S  *  X  )) 

AX  =  -0.25  /  (X  *  Al) 

•  BX  =  0.5  /  SQRT (S  *  X) 

O  *  -T  /  (  A  *  B) 

F  =  -X  ✓  (S  *  A  *  B) 

DX  s  ((AX  *  B  ♦  BX  •  A)  *  X  ♦  A  *  B)  /  (4.0  *  S  *  A  *  A  • 

1  8  *  3  *  X  *  X) 

^  FX  =  -(A  •  B  -  (AX  *  B  ♦  BX  *  A)  *  X)/(S  *A*A*B*B> 

RETURN 
END 
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Table  5 


■»# 


SUBROUTINE  UNORTH 


SUBROUTINE  UNORTH  (PROP,  X,  Y1 
OIMENSION  PR0P( 3) 

PR0P(1)  =  X  COORDINATE  OF  BOUNDARY  LINE  1 

PROP ( 2 )  =  X  COORDINATE  OF  BOUNDARY  LINE  3 

PROP( 3)  =  Y  COORDINATE  OF  BOJNO ARY  LINE  4 

PROP ( 4 1  s  Y  COORDINATE  OF  BOUNDARY  LINE  2 

COMMON  /  FQA  /  NSHELL ,  A,  B,  AX,  AY,  8X,  BY,  XG,  YG,  ZG 

COMMON  /  FQB  /  C,  CX,  CY,  D,  OX,  DY,  E,  EX,  EY,  F,  FX,  FY 

elliptical  cone 

»**»«»«*  '♦*»#♦**♦*♦*♦**•#»•**',♦♦»»♦•*#•******#*****♦**♦*»» 

X  s  W  *  S  *  C0S(T»,  Y  =  V  *  S  *  SIN ( T > ,  Z  =  S 
M  *  PROP (5) ,  SEMI-MAJOR  AXIS  WHEN  S  =  1,  INPUT  ON  CARD  G-2 
V  =  PROP ( 6) ,  SEMI-MINOR  AXIS  WHEN  S  *  1,  INPUT  ON  CARO  G-2 
T  IS  ANGLE  IN  ELLIPTICAL  COORDINATES,  RADIANS 
PROP(  3 )  AND  PROP  (41  ARE  IN  RADIANS  FOR  THIS  EXAMPLE 


X  s  X  ♦  PR0P(1) 

Y1  =  Y  ♦  PROP(3) 

XG  =  W  *  X  *  COS(Yl) 

YG  =  V  *  X  *  SIN(Y1) 

ZG  =  X 

T1  =  ATAN( (PROP (5)  /  PR0P(6>)  *  TAN(Yll) 

M2  s  PROP (5 )  *  PROP ( 5 ) 

V2  =  PROP(6»  *  PR0P( 6) 

MV  =  PROP(5>  *  PROP ( 6 ) 

SH  =  S I N (Til 
CH  =  COS(Tl) 

SH2  =  SH  *  SH 
CH2  =  CH  *  CH 

A  =  SORT  (1.0  H2  *  CH2  ♦  V 2  *  SH2) 

B  =  X  *  SORT ( W2  *  SH2  ♦  V2  *  CH2 ) 

C  =  -X  *  (W2  -  V 2)  *  SH  *  CH 
AY  =  -(M2  -  V2>  *  SH  *  CH  /  A 
BX  =  B  /  X 

BY  =  X  *  X  *  SH  *  CH  *  (M2  -  V 2)  /  3 

CX  =  C  /  X 

CY  =  -X  *  (W2  -  V2 1  *  (CH2  -  SH2 ) 

H  =  X  *  SQRT ( M2  *  SH2  ♦  V2  *  CH2  ♦  M2  *  V2> 

F=MV*X*X/H 

FX  =  F  /  X 

FY  =  -MV  *  X  *  SH  *  CH  *  (H2  -  V2 >  *  (X  /  H) **3 


RETURN 

ENO 
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Table  6 

SUBROUTINE  TEMP 


SUBROUTINE  TEMP  (X  ,  Y , T , API , AP2> 


*0 


.T  *  0  •  0 

50.»(l.-.b»COS<Y*J. > 
API*. 00001 
AP2=.30001 
'"RETURH 
END 


,19/UO.)  ) 


Table  7 

SUBROUTINE  WALL 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


SUBROUTINE  HALL  (X,  Y,  CCC) 

OIMENSION  CCC (6(6l 
COMMON  /  OFST  /  TO,  l 

l  -  OISTANCE  FROM  REFERENCE  SURFACE  TO  MIP'URFACS  OF  SHELL 

HALL 

TO  =  TOTAL  thickness  of  shell  hall 


C  F  8  2 

COMMON  /  MONO  /  AT,  EX1,  XN1,  EYi,  G 


C  F  6  3 

COMMON  /  SKEM  /  E3,  U3,  T 3 ,  TM,  A,  8,  M3,  AK3 


C  F  8  4 

COMMON  /F I8R/  EF.EM, UF,  UN, LAYERS, TT  C  20  » ,XX(20» , BE  120) , 01 ?0 ) 


C  F  8  5 

COMMON  /LAYOl/  TL  ( 2:  (  ,  E X  5  <  20  I  , £Y  5  ( 20 )  , UX Y  «  20 )  ,  G5 (20  >  , L AY S 


C  F  8  6 

COMMON  /  CORR  /  CTb,  £6,  Ub,  CCb,  CMb,  COb,  C86 


C  F  6  7 

COMMON  /CORS/  CT7,E7,U7,CC7,CH7,CD7,CB7,ES,US,TS,PHI,ANC 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

C 


OPTION  A 

GENERATES  A  b»b  MATRIX  CCC  ACCORDING  TO  SHELL  MALL 
CONSTRUCTION.  SET  2  AND  COMMON  VARIABLES  FOR  CFBN.  THEN 
CALL  CFB(N) . 


OPTION  3 

SET  Z,  TO,  ANO  THE  6*6  MATRIX  CCC. 
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Table  7  (Cont. ) 
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1 


o 


Section  6 

INPUT  DESCRIPTION 


There  are  two  main  branches  of  the  program,  one  with  and  one  without  plasticity. 

The  plasticity  branch  cannot  be  used  with  variable  material  properties  and  is  restricted 
to  isotropic  wall  construction.  Figure  G-l  shows  the  data  deck  format.  Table  8  is  a 
minimanual  summarizing  the  input  cards.  Figure  6-2  shows  the  sign  convention  for 
stress  and  moment  resultants. 


Fig.  6-1  STAGS  Program  -Data  Deck  Format 


G-l 


Tabic  8 


STAGS  PROGRAM  MINIMANUAL 


INPUT  DESCRIPTION 


Item 

Symbol 

Format 

C-l 

COMENT  (I),  1=1,12 

12A6 

G-l 

NSHELL,  INDIC ,  NLOAD,  NCIIK 

415 

G-2 

PROP  (I),  1=1,8 

8E10.6 

D-l 

NR,  NO ,  NRW1 ,  NRW2 ,  NCL1,  NSTFF 

(INCLUDE  ITEMS  D-2  TO  D-4  ONLY  IF  NR.  EQ.O) 

615 

D-2 

NNX 

15 

D-3 

SEGLX(I),  1=1,  NNX 

8E10.6 

D-4 

NSEGX(I) ,  1=1,  NNX 

(INCLUDE  ITEMS  D-5  TO  D-7  ONLY  IF  NC.  EQ.O) 

1615 

D-5 

NNY 

15 

D-6 

SEGLY(J),  J=l, NNY 

8E10.6 

D-7 

NSEGY(J),  J=l,  NNY 

(INCLUDE  ITEM  D-8  ONLY  IF  NR.LT.O) 

1615 

D-8 

X(I),  1=1,  NR 

(INCLUDE  ITEM  D-9  ONLY  IF  NC.  LT.O) 

8E10.6 

D-9 

Y(J),  J=1,NC 

(INCLUDE  ITEMS  1-1  TO  1-4  ONLY  IF  INDIC.  EQ.3) 

8E10.6 

1-1 

AE ,  XNU ,  AT ,  AK2 

4  E10.6 

1-2 

NL.IC 

215 

1-3 

S(I),  1=1,  ic 

8E10.6 

1-4 

E(I),  1=1,  IC 

8E10.6 

B-l 

IBLN(I) ,  1=1,4 

(INCLUDE  ITEM  B-2  ONLY  IF  IBLN(I).  EQ.  0) 

415 

B-2 

ICOND(I) ,  1=1,4 

415 

L-l 

NN(K),  LFLG(K),STLD(K),LSTP(K),  MXL(K) 

2I5.3E10.6 

L-2 

PZ,PY,PX,JZ,  JY,  JX,  L,M 

(REPEAT  ITEM  L-2  NN  TIMES) 

(REPEAT  ITEMS  L-l  AND  L-2  FOR  K=l,2) 

(IF  NO  DISCRETE  STIFFENERS  ARE  PRESENT,  GO 

TO  P-1A  CARD  FOR  BIFURCATION,  P-1B  CARD  FOR 
NONLINEAR,  AND  0-1  CARD  FOR  LINEAR  STRESS 
ANALYSIS) 

3E10.6.5I5 

S-I 

IRGS,  ITRN,  IRSO,  ISTR,  ITSN,  ISSO 

(INCLUDE  ITEMS  S-2  TO  S-4  ONLY  IF  IRGS.GT.O) 

615 

S-2 

IRN(I) ,  IRTP(I) ,  IRNA(I) ,  IRNB(I) , XRN(I) ,  Y 1  RN(I) ,  Y2RN(I) 
(REPEAT  ITEM  S-2  FORI  1,  IRGS) 

415, 3E 10. 6 

S-3 

ERN(J) ,  ZARN(J) ,  ZIXRN(J) ,  ZIZRN(J) ,  ZJRN(J) ,  EZRN(J) ,  ZK1 
(REPEAT  ITEM  S-3  FOR  J-l, ITRN) 

(INCLUDE  ITEM  S-4  ONLY  IF  IRSO.GT.O) 

7E10.6 

S-4 

Z1,X1,Z2,X2,  Z3,X3,Z4,X4 

(REPEAT  ITEM  S-4  FOR  J- 1 ,  ITRN) 

(INCLUDE  ITEMS  S-5  TO  S-7  ONLY  IF  ISTR.GT.O) 

8E10.6 

6-2 


Table  8  (Cont.) 

Item  Symbol  Format 


S-5 

ISN(I) ,  ISTP(I) ,  ISTA(I) ,  ISTB(I) ,  YSN(I) ,  X1SN(I)  ,X2SN(I) 
(REPEAT  ITEM  S-5  FOR  I=1,ISTR) 

4I5.3E10.6 

S-6 

ESN(J)  ,XASN(J) ,  XIYSN(J) ,  XIZSN(J) ,  XJSN(J) ,  E ZSN(J) ,  XK1 
(REPEAT  ITEM  S-6  FOR  J=1,ITSN) 

(INCLUDE  ITEM  S-7  ONLY  IF  ISSO.GT.O) 

7E10.6 

S-7 

Z1,Y1,Z2,Y2,Z3,Y3,Z4,Y4 

(REPEAT  ITEM  S-7  FOR  J=1,ITSN) 

(INCLUDE  ITEMS  P-1A  TO  P-1A2  ONLY  FOR 
BIFURCATION  ANALYSIS) 

8E10.6 

P-1A 

DE  LBIF ,  SHIFT ,  IBOND ,  ISHIFT ,  ITERAT 

(INCLUDE  ITEM  P-1A1  ONLY  IF  IBOND.  EQ.  1) 

2E10.6.3I5 

P-1A1 

JBLN(I),  1=1,4 

415 

P-1A2 

JCOND(I),  1=1,4 

(INCLUDE  ITEM  P-1B  ONLY  FOR  NONLINEAR 

ANALYSIS) 

415 

P-1B 

DELX,  WUND,  ISTART,  ISEC ,  ICUT,  INEWT,  ISTRAT 

2E10.6.5I5 

0-1 

IPR1,  IPR2 ,  IPR3 ,  IPX ,  IPY ,  IPRD,  IPRS,  IP  LOT,  IHIST 
(INCLUDE  ITEM  0-2  ONLY  IF  IPX.GT.O) 

915 

0-2 

XPL(I),  1=1,  IPX 

(INCLUDE  ITEM  0-3  ONLY  IF  IPY.GT.O) 

8E10.6 

0-3 

YPL(I),  1=1,  IPY 

(INCLUDE  ITEM  0-4  ONLY  IF  IHIST.  GT.O) 

8E10.6 

0-4 

XI,  Y1.X2.Y2.X3,  Y3.X4,  Y4 

(IF  INDIC.EQ.3  NO  MORE  CARDS  ARE  REQUIRED) 

8E10.6 

M-l 

IWALL,  NSTRI,  NRING,  IP,  IM,  JM 

(INCLUDE  ITEM  M-2B  ONLY  IF  IWALL.  EQ.  2) 

615 

M-2B 

AT,  EX1.XNU,  Z,  EY1,  G 

(INCLUDE  ITEMS  M-2C1  AND  M-2C2  ONLY  IF 

IWALL.  EQ.  3) 

6E10.6 

M-2C1 

T3,E3,U3,TH,  A,  B,H3,  AK3 

8E10.6 

M-2C2 

Z 

(INCLUDE  ITEMS  M-2D1  AND  M-2D2  ONLY  IF 

IWALL.  EQ.  4) 

E10.6 

M-2D1 

EF,  EM,  UF,  UM,  Z,  LAYERS 

5E10.6.I5 

M-2D2 

TT(J),XX(J),BE(J),  O(J) 

(REPEAT  ITEM  M-2D2  FOR  J=l,  LAYERS) 

(INCLUDE  ITEMS  M-2E1  AND  M-2E2  ONLY  IF 

IWALL.  EQ.  5) 

4E10.6 

M-2E1 

Z, LAYS 

E10.6.I5 

M-2E2 

TL(J) ,  EX5(J) ,  EY5(J) ,  UXY(J) ,  G5(J) 

(REPEAT  ITEM  M-2E2  FOR  J=l,  LAYS) 

(INCLUDE  ITEM  M-2F  ONLY  IF  IWALL.  EQ.  6) 

5E10.6 

M-2F 

CT6,  E6,  UG,  CC6,  CH6,  CD6,  CB6,  Z 

(INCLUDE  ITEMS  M-2G1  AND  M-2G2  ONLY  IF 
IWvLL.EQ.7) 

8E10.6 

M-2G1 

CT7,  i^7,  U7,  CC7,  CH7.  CD7,  CB7  Z 

8E10.  6 

6-3 


Item 

M-2G2 

N-l 

N-2A 

N-2B 

N-3 

N-4A 

N-4B 

V-l 

V-2 

V-3 

V-4 


Table  8  (Cont.) 


Symbol 


TS,  ES,  US,  PHI,  ANC 

(INCLUDE  ITEMS  N-l  TO  N-2B  ONLY  IF  NSTRI.EQ.  1) 
El,  U1 ,  Oil ,  Dl,  AK1 

(INCLUDE  ITEM  N-2A  ONLY  IF  OI1.EQ.O) 

Tl.Hl 

(INCLUDE  ITEM  N-2B  ONLY  IF  OI1.EQ.  1) 

Al,  SI1.XI1,  SI,  EZ1,  HI 

(INCLUDE  ITEMS  N-3  TO  N-4B  ONLY  IF  NRING.EQ.l) 
E2,U2,  012,  D2,  AK2 

(INCLUDE  ITEM  N-4A  ONLY  IF  OI2.EQ.O) 

T2,H2 

(INCLUDE  ITEM  N-4B  ONLY  IF  012.  EQ.  1) 

A2,  SI2.XI2,  S2,  EZ2.H2 

(INCLUDE  ITEMS  V-l  TO  V-4  ONLY  IF  IWALL.EQ.9) 
TD,  Z 

XM(I),  1=1,  IM 
YM(J) ,  J=l, JM 

((TDEG(L, M, N),  EX(L, M,  N),  EY(L, M,  N),  U(L,  M,N), 

G(L,  M,  N,),  A1(L, M, N),  A2(L,  M,  N),  L=l, IP) ,  M-l.MSTA) 
WHERE  MSTA=IM  IF  NC.GT.NR,  ELSE  MSTA=JM. 
(REPEAT  ITEM  V-4  FOR  N=1,NSTA,  WHERE 
NSTA=JM  IF  NC.GT.NR,  ELSE  NSTA=IM) 


Format 

5E10.6 

5E10.6 

2E10.6 

6E10.6 

5E10.6 

2E10.6 

6E10.6 

2E10.6 

8E10.6 

8E10.6 

7E10.6 
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STRESS  RESULTANTS 


T  "vx 


MOMENT  RESULTANTS 


Fig.  6-2  Sign  Convention  for  Stress  and  Moment  Resultants 


C-l  Comment  Card 

The  Comment  or  Case  Title  card  may  contain  any  Hollerith  Text.  This  comment  is 
printed  at  the  beginning  of  the  output  for  the  case. 

Variable  Format  Columns  Description 

COMENT  12A6  1-72  Case  Title 


GEOMETRY 


-1  Shell  and  Analysis  Type  Definition  Card 


The  card  is  used  to  define  the  surface  type  of  the  shell  by  a  single  intege:  •  itry.  It 
also  serves  to  indicate  the  type  of  analysis  desired. 


Variable 

Format 

Columns 

Description 

NSHELL 

15 

1-5 

Shell  type  may  vary  from  1  through  N  as 
described  under  G-2  card. 

NSHELL  =  1,  Cylinder 

2,  Cone/ Annular  Plate 

3,  Plate 

4,  Sphere 

5,  Paraboloid 

6,  Elliptic  Cylinder 

7,  Ellipsoid 

8,  Torus 

9,  Hyperboloid 
10,  Elliptic  Cone 

11, 12,  User  written  subroutine 

0  -  Linear  solution  only 

1  -  Bifurcation  analysis 

2  -  Nonlinear  analysis  (elastic) 

3  -  Nonlinear  analysis  (plastic) 

Number  of  load  systems  to  be  applied 
independently  (presently  NLOAD  is 
restricted  to  1  or  2)  (Section  3.7, 
Loading;.  Note:  Load  cards  for  system  B 
omitted  if  NLOAD  *  1. 

NCHK  15  16-20  0  -  Execute  (no  3D  plot) 

1  -  Do  not  execute.  Provide  input  data 

check  and  3D  plot 

2  -  Execute  and  provide  3D  plot 


INDIC  15  6-10 


NLOAD  15  11-15 


6-6 


G-2  Surface  Constants  Card 


This  card  will  contain  the  various  measurements  related  to  the  particular  type  of  sur¬ 
face  defined  by  the  NSHELL  integer  on  the  G-l  card. 


Variable  Format  Columns  Description 

PROP  (I)  8E  10.6  1-80  Surface  properties  depending  on  NSHELL 

(Fig.  6-3).  All  angles  in  degrees  for 
NSHELL  less  than  11.  See  Table  4, 
p.  5-7,  for  definition  of  PROP  (I)  when 
NSHELL  greater  than  10. 


Note:  X-|  and  Y-i?  coordinates  designation  are  synonymous  for  all  shell  types. 


CYLINDER 

|W3HCU  •  I  I  f*.*1 


CONE  ANNULAR  PLATE 

[nSHEU ■  2  1  R*  R* 


Fig.  6-3  Types  of  Shell  Surfaces  Defined  by  NSHELL  Integer 


PLATE 


NSHELL= 3 


Note:  x-y  and 

X-Y  are  not  synonymous. 

(See  Nomenclature.) 


NSHELL *  5 


Prop( I ) «  L 
Prop(2)=  6  (deg) 
Prop(3)s  R 
Prop{4)=  R, 


PARABLOID 
y2*  zz  *  4Rx 


Fig.  6-3  (Cont.) 
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Prop(  I )  *  L 
Prop(2)  *  0L 
Prop(3)  = 
Prop(4)*  Rz 
Prop(5)  *  Ry 


Prop(l)*CL  (defl ) 
Prop(2)»fu  (deg) 
Prop(3)*0  (deg) 
Prop(4)«  Ry 
Prop(5)  ■  Rx 


Prop(  I ) «  £L  (deg) 
Prop(2)*£y  (deg) 
Prop(3)c  8  (deg) 
Prop(4)  *  R 
Prop(5)*  R, 


Fig.  6-3  (Cont.) 
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HYPERBOLOID 


Prop(l)*£L  (deg) 
Prop(2)*fu  (deg) 
Prop(3)*6  (deg) 
Prop(4)-  R 
Prop(5)*  R| 


x  (Axi*  of  Revolution) 


Prop{  I )  •  L 
Prop<2)  •  L. 
Prop(3)  •  6l  (d«*) 
Prop(4)  •  Bu  (dtg) 
Prop!  3)  ■  R, 
Prop<6)  ■  Ry 


ELLIPTIC  CONE 


I 


Fig.  6-3  (Cont.) 

NSHELL  =  11 

User  supplied  subroutine  (ORTH)  for  shells  described  by  orthogonal  surface  coordinate 
lines  and  not  included  in  NSHELL  =  1  through  10  (see  Section  5.4). 


NSHELL  =  12 

User  supplied  subroutine  (UNORTH)  for  shells  described  by  nonorthogonal  surface 
coordinate  lines  and  not  included  in  NSHELL  =  1  through  10  (see  Section  5.5). 
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DISCRETIZATION 


D-l  Mesh  Definition  Card 


Variable 

Format 

Columns 

Description 

NR 

15 

1-5 

Number  of  rows  (Stations  along  the  X 
coordinate).  If  NR  negative,  spacing  is 
variable,  read  in  X  coordinates  on 

Card  D-8.  If  NR  =  0,  spacing  is  con¬ 
stant  within  each  of  a  number  of  segments 
but  varies  from  one  segment  to  another 
(see  D-2  card). 

NC 

15 

6-10 

Number  of  columns  (Stations  along  the  V 
coordinate'.  If  NC  negative,  spacing  is 
variable,  read  in  Y  coordinates  on 

Card  D-9.  If  NC  =  0,  spacing  is  con¬ 
stant  within  each  of  a  number  of  segments 
but  varies  from  one  segment  to  another. 
(See  D-5  card.) 

NRW1 

15 

11-15 

Row  number  of  one  edge  of  cutout. 

NRW2 

15 

16-20 

Row  number  of  other  edge  of  cutout. 

NCL1 

15 

21-25 

Column  number  of  edge  of  cutout. 

NSTFF 

15 

26-30 

0  -  No  discrete  stiffeners 

1  -  Discrete  stiffeners  (input  data  on 

S-l  card) 

(This  card  continued  next  page.) 
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Note: 


NRW1,  NRW2,  and  NCL1  define  a  rectangular  cutout  adjacent  to  boundary 
line  4.  In  the  case  of  no  cutout,  NRW1,  NRW2,  and  NCL1  should  be  blank. 
Additional  cutouts  or  cutouts  of  more  general  shape  can  be  included  through 
specification  of  a  zero  modulus  of  elasticity  in  the  appropriate  area  (use 
IWALL  =  1  on  M-l  card). 


If  NR  and  NC  are  positive  numbers,  omit  cards  D-2  through  D-9. 
D-2  X-Segment  Card 


This  card  should  be  included  only  if  NR  is  zero  on  the  D-l  card. 

Variable  Format  Columns  Description 

NNX  15  1~5  Number  of  segments  in  X  direction  with 

constant  spacing. 


D-3  X-Segment  Length  Definition  Cards 


These  cards  should  be  included  only  if  NR  is  zero  on  the  D-l  card. 


Variable 

SEGLX(I) 

I  =  1,  NNX 


Format  Columns 


0E  10.6 


1-80 


Description 

"Length"  of  Segment  I,  (i.e.,  difference 
between  extreme  values  of  X  coordinate) . 


D-4  X-Segment  Spacing  Definition  Cards 


These  cards  should  be  included  only  if  NR  is  zero  on  the  D-l  card. 


Variable 

NSEGX(I) 

I  =  1,  NNX 


Format  Columns 


16  15 


1-80 


Description 

Number  of  mesh  spaces  within  Segment  I. 
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D-5  Y-Segment  Card 


This  card  should  be  included  only  if  NC  is  zero  on  the  D-l  card. 


Variable 

Format 

Columns 

Description 

NNY 

15 

1-5 

Number  of  segments  in  Y  direction  with 
constant  spacing. 

D-6  Y-Segment  Length  Definition  C. 

These  cards  should  be  included  only  if  NC  is  zero  on  the  D-l  card. 

Variable  Format  Columns  Description 

SEGLY(J)  8E  10.6  1-80  ’’Length"  of  Segment  J  (i.  e. ,  difference 

J  =  lrNNY  between  extreme  values  of  Y  coordinate). 


D-7  Y-Segment  Spacing  Definition  Cards 

These  cards  should  be  included  only  if  NC  is  zero  on  the  D-l  card. 


Variable  Format  Columns 

NSEGY(J)  16  15  1-80 

J  =  l.NNY 

Unless  NR  or  NC  is  negative,  omit  D-8  and  D-9  cards. 


Description 

Number  of  mesh  spaces  within  Segment  J. 


D-8  X-Coordinate  Cards 


These  cards  should  be  included  only  if  NR  is  negative  on  the  D-l  card. 


Variable  Format  Columns  Description 

X(I)  8E  10.6  1-80  X  coordinate  for  Row  I  (Must  be 

(I  =  1,NR)  monotonically  increasing) 
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D-9  Y -Coordinate  Cards 


These  cards  should  be  included  only  if  NC  is  negative  on  the  D-l  card. 


Variable  Format  Columns  Description 

Y(J)  8E  10.6  1-80  Y  coordinate  for  Column  J  (Must  be 

(J  =  1,NC)  monotonically  increasing) 

Note: 

Unless  the  plasticity  branch  is  used,  omit  cards  1-1  through  1-4,  and  go  to  the 
B-l  card. 
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PLASTICITY 


Cards  1-1  through  1-4  are  used  only  for  the  plasticity  branch,  i.e. ,  INDIC  =  3  on  the 
G-l  card.  For  the  elastic  branch  go  to  the  B-l  card. 

1-1  Plasticity  Definition  Card-1 


Variable 

Forma* 

(  olumns 

Description 

AE 

E  10.6 

1-10 

Young's  Modulus 

XNU 

E  10.6 

11-20 

Poisson's  Ratio 

AT 

E  10.6 

21-30 

Shell  Thickness 

AK2 

E  10.6 

31-40 

Square  of  the  ellipse-ratio  of  yield 
surface.  Usually  AK2  =  3.0 

If  the  elastic  branch  is  used, corresponding  information  is  read  at  a  different  place. 


1-2  Plasticity  Definition  Card -2 


Variable 

Format 

C  .umns 

NL 

15 

1-5 

Number  of  points  across  wall  thickness. 
Must  be  odd  number  and  not  less  than  3 
or  more  than  9. 

IC 

15 

6-10 

Number  of  material  components  (Number 
of  points  defined  on  the  stress- strain 

• 

curve).  (Fig.  6-4)  IC  <  10 

Fig.  6-4  Stress-Strain  Curve 
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Variable  Format,  Columns, 


S(I) 

I  --  1,  IC 


E  10.  6 


1-4  Plasticity  Definition  Card-4 
Variable  Formal  Cohmm§- 


E(I) 

I  =  1,IC 


E  10.  6 


Description 

Stress  values  on  the  stress-strain  curve. 
See  Fig.  6-4. 


Description 


Corresponding  strain  values ^  on  the 
stress-strain  curve.  See  Fig.  6-4. 


Note,  E(l)  need  not  be  punched,  It  will  be  computed  by  the  program  as 

E(l)  =  S(1)/AE 
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BOUNDARY  CONDITION 


B-l  Boundary  Condition  Card-1 


Do  not  specify  loads  and  displacements  which  are  in  conflict,  such  as  simple  support 
and  specified  tangential  displacement. 


For  bifurcation  buckling  (INDIC  =  1),  the  conditions  defined  here  will  apply  to  the  pre- 
buckling  displacements.  If  the  boundary  conditions  for  incremental  displacement  are 
different,  they  will  be  specified  on  the  P-1A,  P-1A1,  and  P-1A2  cards  below. 


Boundary  lines  are  numbered  1  through  4  according  to  the  following  convention: 


Row  1  -  Boundary  line  1 

Column  NC  -  Boundary  line  2 


Row  NR 
Column  1 


-  Boundary  line  3 

-  Boundary  line  4 


Note:  Side  1  X  =  0,  Side  4 


Variable 

1BLN(I) 

I  =  1.4 


Format 


4  15 


Y  =  0. 
Columns 

1-20 


Y  or  0,  r) 


Description 

Boundary  code  for  Line  I. 

0  “  Specified  on  B-2  cards 

1  -  Simple  support  (defined  below) 

2  -  Clamped  (defined  below) 

3  -  Unrestrained 

4  -  Symmetry  (defined  below) 

5  -  Anti-symmetry  (defined  below) 

6-  Closed  shell  in  the  rj  direction  (see 
comment  below  and  pages  2-2  and 
2-3). 


Here  simple  support  means  w  =  v 

w  =  u 

(This  card  continued  next  page. ) 


0  on  lines  1  and  3 
0  on  lines  2  and  4 
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clamped  means  that  all  displacements  and  rotations  are  restrained  and  cannot  be 

used  for  loaded  edge 

Symmetry  implies  u 

v 

Antisymmetry  implies  v 

u 

The  option  IBLN  =  6  indicates  a  closed  shell.  In  that  case,  set  IBLN(2)  =  IBLN(4)  =  6. 
Do  not  set  IBLN  (1)  and  IBLN(3)  equal  to  six.  Boundary  lines  of  aero  length  (e.g. , 
apex)  have  not  been  incorporated  in  the  program. 


=  0  =  0  on  lines  1  and  3 

=  0  =  0  on  lines  2  and  4 

=  w  =  u,  =  0  on  lines  1  and  3 

=  w  v, x  =  0  on  lines  2  and  4 

*  ir 


B-2  Boundary  Condition  Cards.  Use  only  if  IBLN(I)  =  0  on  the  B-l  card.  One  B-2 
card  for  each  I  with  IBLN  =  0.  The  cards  refer  to  boundary  lines  with  increasing 
numbers. 


Variable  Format  Columns 

ICOND(J)  4  15  1-20 

J  =  1,  4 


Description 

Freedom  of  movement  of  Line  I  in  regard 
to  the  w  ,  v,  u  ,  and  0  displacements, 
where  w  ,  v  ,  and  u  are  displacements  in 
the  direction  of  the  coordinates  Z,  Y,  and 
X,  respectively,  and  0  is  the  rotation 
around  a  tangent  to  the  edge 

1  -  Free  to  move 

0  -  No  movement 


Y.tj.v.O 


X,{  ,u 
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LOADS 


A  Load  System  is  defined  by  an  L-l  card  followed  by  the  appropriate  L-2  cards.  If 
NLOAD  on  the  G-l  card  is  2,  two  Load  Systems  are  required  -  one  for  Load  System  A 
and  another  for  Load  System  B.  The  input  data  cards  for  Load  System  A  must  be  com¬ 
pleted  before  the  input  data  cards  for  Load  System  B  begin  .  The  use  of  Load  System 
B  in  bifurcation  analysis  is  explained  in  Section  3.7,  Loading. 

A  base  load  for  each  Load  System  is  defined  by  the  set  of  L-2  cards  together  with  any 
loads  that  may  be  introduced  by  a  user-wrtten  subroutine  (USRLD).  STLD,  LSTP, 
and  MXL  refer  to  the  base  load  of  the  system. 


L-l  Card 


Variable 

Format 

Columns 

Description 

NN 

15 

1  -5 

Number  of  L-2  cards  required  to  describe 
the  Load  System. 

LFLG 

15 

6-10 

User -Load  Flag. 

0  -  User  does  not  have  own  subroutine  to 
define  Load  System. 

1  -  User  has  own  subroutine  to  define 

Load  System. 

STLD 

E  10.6 

11  -20 

Starting  Load  Factor.  (The  initial  load  is 
STLD  times  the  base  load.)  For  nonlinear 
analysis,  STLD  is  the  current  starting 
Load  Factor.  (See  sample  case  1  Second 
Run,  Page  7-10).  For  linear  analysis, 
the  total  load  equals  the  initial  load.  For 
bifurcation  analysis,  the  critical  load  is 
the  eignevalue  times  the  initial  load. 

LSTP 

E  10.6 

21  -30 

Load  Step  Increment.  (The  load  incre¬ 
ment  is  LSTP  times  the  base  load.) 
Meaningful  only  for  nonlinear  analysis. 

MXL 

E  10.6 

31  -40 

Maximum  Load.  (The  maximum  load  is 
MXL  times  the  base  load.)  Meaningful 
only  for  nonlinear  analysis.  May  be  used 
to  freeze  either  Load  System  A  or  Load 
System  B. 
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L-2  Cards 


Do  not  prescribe  loads  and  displacements  on  the  same  L-2  card. 


Variable 

Format 

Columns 

PZ.PY.PX 

3  E  10.6 

1-30 

JZ, JY.JX 

3  15 

31  -45 

L,  M  2  15  46  -55 


Description 

Load  Element.  (Positive  outward  and 
toward  increasing  X  and  Y.) 

Refer  to  PZ,PY,  and  PX,  respectively 
-1  -  Displacement 
0  -  Omit 

1  -  Point  Force 

2  -  Line  Load  along  Row 

3  -  Line  Load  along  Column 

4  -  Pressure  Load 

5  -  Live  Pressure  Load 

(Use  only  for  L  =  M  =  0. ) 

Row  and  Column  number  of  the  mesh- 
point  where  the  particular  load  element 
is  applied. 


Note: 

When  the  row  number  is  entered  as  zero,  the  load  element  is  assumed  to  act 
at  every  mesh-point  on  the  column  indicated  by  the  column  number.  If  the 
column  number  is  entered  as  zero,  the  load  Element  is  assumed  to  act  at 
every  mesh -point  on  the  row  indicated  by  the  row  number.  If  both  row  and 
column  numbers  are  zero,  then  the  load  element  is  assumed  to  act  at  each 
mesh-point. 


For  inplane  displacements,  only  uniform  displacement  of  one  of  the  four 
boundary  lines  may  be  prescribed.  That  is,  either  L  or  M  must  be  zero. 


The  energy  method  on  which  the  program  is  based  assumes  that  loads  and  structure 
comprise  a  conservative  system.  Because  there  is  some  controversy  about  the  re¬ 
quirements  for  a  system  with  variable  load  to  be  conservative,  the  live  load  option  is 
internally  suppressed  unless  the  pressure  is  uniform. 

Here  follows  L-l  and  L-2  cards  for  Load  System  B  if  NLOAD  (G-l  card)  is  equal  to  2. 
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DISCRETE  STIFFENERS 


S-l  Discrete  Stiffener  Definition  Card  {Ring  is  defined  as  a  stiffener  on  constant 
X-coordinate;  stringer  is  defined  as  a  stiffener  on  constant  Y -coordinate. )  Note  that 
"smeared"  stiffener  can  be  defined  in  lieu  of  or  in  addition  to  the  discrete  stiffeners 
(N-l  through  N-4B  cards  below). 

Note: 


If  no  discrete  stiffeners  are  present,  NSTFF  equal  zero  on  D-l  card;  go  to  the 
P-1A  card  for  bifurcation  analysis,  to  P-1B  card  for  nonlinear  analysis,  and 
to  the  0-1  card  for  linear  stress  analysis. 


Variable 

Format 

CoUimns 

Description 

IRGS 

15 

1-5 

Number  of  rings,  IRGS  <  80 

ITRN 

15 

6-10 

Number  of  distinct  type  of  rings. 

ITRN  £  30 

mso 

15 

11-15 

Stress  output  will  be  given  at  IRSO  points 
In  the  ring  cross-section.  IRSO  <  4 

ISTR 

15 

16-20 

Number  of  stringers.  ISTR  £  80 

itsn 

15 

21-25 

Number  of  distinct  type  of  stringers. 

ITSN  <  30 

ISSO 

15 

26-30 

Stress  output  will  be  given  at  ISSO  points 
in  the  stringer  cross-section.  ISSO  <  4 

If  there  are  no  rings  (IRGS  =  0),  go  to  the  S-5  card. 
S-2  Ring  Delimiter  Cards 


Include  these  cards  only  if  IRGS  on  Card  S-l  is  not  0. 


Variable 

Format 

Columns 

Description 

IRN(I) 

15 

1-5 

Row  number  of  ring 

IRTP(I) 

15 

6-10 

Ring  type 

IRNA(I) 

15 

11-15 

Starting  column  of  ring 

(This  card  continued  next  page.) 

f 

I 
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Variable 

Format 

Columns 

Description 

ERNB(I) 

15 

16-20 

Ending  column  of  ring.  (If  IRN,  IRNA, 
and  IRNB  are  set  equal  to  zero,  the  posi 
tion  of  the  ring  is  given  by  coordinate 
values  rather  than  by  row  and  column 
numbers) 

XRN(I) 

E  10.6 

21-30 

X  coordinate  of  ring 

YIRNfl) 

E  10.6 

31-40 

Y  coordinate  of  start  of  ring 

Y2RN(I) 

E  10.6 

41-50 

Y  coordinate  of  end  of  ring 

Notes: 

1.  XRN,  Y1RN,  and  Y2RN  are  optional  data  elements  and  will  be  used  to 
define  ring  location  only  if  IRN,  ERNA,  and  IRNB  are  zero,  respectively. 

2.  If  XRN,  Y1RN,  and  Y2RN  do  not  coincide  with  a  grid  line,  the  ring  will 
be  placed  (by  the  program)  at  the  closest  grid  line. 

3.  Repeat  this  card  for  I  =  1,  ERGS. 

S-3  Ring  Description  Cards 


Include  these  cards  only  if  IRGS  card  S-l  is  not  zero.  Card  format  is  8E  10.6.  One 
card  is  required  for  each  ring  type  (ITRN  cards). 


(This  card  continued  next  page.) 
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Variable 

Format 

Columns 

ERN(J) 

E  10.6 

1-10 

ZARN(J) 

E  10.6 

11-20 

ZIXRN(J) 

E  10.6 

21-30 

ZIZRN(J) 

E  10.6 

31-40 

ZJRN(J) 

E  10.6 

41-50 

EZRN(J) 

E  10.6 

51-60 

ZK1 

E  10.6 

61  -  70 

Note: 

Repeat  above  data  for  J  =  1,  ITRN. 


Description 

Young's  Modulus  for  ring  type  J. 

Cross-Section  Area  for  ring  type  J. 

Moment  of  Inertia  about  X'-axis  for  ring 
type  J. 

Moment  of  Inertia  about  Z'-axis  for  ring 
type  J. 

Torsional  stiffness  GJ  for  ring  type  J 

Eccentricity  in  Z'  direction  for  ring 
type  J.  Outside  ring  —  distance  from 
outer  shell  surface  to  ring  centroid 
(positive  in  +  Z  direction).  Inside 
ring  -  distance  from  inner  shell 
surface  to  centroid  (positive  in  -  Z 
direction) . 

0.  -  internal  rings 
1.  -  external  rings 


S-4  Ring  Stress  Output  Card 


Include  these  cards  only  if  IRSO  on  Card  S-l  is  not  0. 


Variable 

Format 

Column 

Description 

Z1 

E  10.6 

1-10 

Z-coordinate  for  the  first  point  with 
stress  output. 

XI 

E  10.6 

11-20 

X'-coordinate  for  the  same  point. 

Z2 

E  10.6 

21-3o"> 

X2 

E  10.6 

31-40 

Same  information  for  other  points  with 

Z3 

E  10.6 

41-50 

>  stress  output.  Read  as  many  pairs  of 

coordinates  as  indicated  by  IRSO  on  the 

X3 

E  10.6 

51-60 

S-l  card. 

Z4 

E  10.6 

61-70 

X4 

E  10.6 

71-80  > 

Note: 

Repeat  above  data  for  each  ring  type,  J  =  1,  ITRN. 
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S-5  Stringer  Delimiter  Cards 

If  there  are  no  stringers  (ISTR  =0),  go  to  the  P-lA  card  for  bifurcation  analysis,  P-1B 
for  nonlinear  analysis  and  to  the  0-1  card  for  linear  stress  analysis. 


Include  these  cards  only  if  ISTR  on  card  S-l  is  not  0. 


Variable 

Format 

Columns 

Description 

ISN  (I) 

15 

1-5 

Column  number  of  stringer 

ISTP(I) 

15 

6-10 

Stringer  type 

ISTA  (I) 

15 

11-15 

Starting  row  of  stringer 

ISTB(I) 

15 

16-20 

Ending  row  of  stringer.  (If  ISN,  TC^'.’A,  and 
ISTB  are  set  equal  to  zero,  the  position  of 
the  stringer  is  given  by  coordinate  values 
rather  than  by  column  and  row  numbers). 

YSN(I) 

E  10.6 

21-30 

Y  coordinate  of  stringer 

X1SN(I) 

E  10.6 

31-40 

X  coordinate  of  start  of  stringer 

X2SN(I) 

E  10.6 

41-50 

X  coordinate  of  end  of  stringer 

NOTES:  1.  YSN,  X1SN,  and  X2SN  are  optional  data  elements  and  will  be  used  tc 
define  stringer  location  only  if  ISN,  ISTA,  and  ISTB  are  zero, 
respectively . 

2.  If  YSN,  X1SN,  and  X2SN  do  not  coincide  with  a  grid  line,  the  strings 
will  be  placed  (by  the  program)  at  the  closest  grid  line. 

3 .  Repeat  this  card  for  1  =  1,  ISTR. 
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S-6  Stringer  Description  Cards 


Include  these  cards  only  if  ISTR  on  card  S-l  is  not  zero.  Card  format  is  8E10.  6. 


Variable 

Format 

Columns 

Description 

ESN  (J) 

E  10.6 

1-10 

Young's  Modulus  for  stringer  type  J. 

XASN(J) 

E  10.6 

11-20 

Cross-section  Area  for  stringer  type  J. 

XIYSN(J) 

E  10.6 

21-30 

Moment  of  Inertia  about  Y*-axis  for 
stringer  type  J. 

XIZSN(J) 

E  10.6 

31-40 

Moment  of  Inertia  about  Z,-axis  for 
stringer  type  J. 

XJSN(J) 

E  10.6 

41-50 

Torsional  stiftness  GJ  for  stringer  type  J 

EZSN(J) 

E  10.6 

51-60 

Eccentricity  in  Z'  direction  for  stringer 
type  J.  Outside  stringer  —  distance  from 
outer  shell  surface  to  stringer  centroid 
(positive  in  +Z  direction).  Inside  stringer  - 
distance  from  inner  shell  surface  to  cen¬ 
troid  (positive  in  -Z  direction). 

XK1 

E  10.6 

61-70 

0.  -  internal  stringers 

1.  -  external  stringers 

Note: 

Repeat 

above  data  for  J 

=  1,  ITSN. 

S-7  Stringer  Stress  Output  Card 

Include  these  cards  only  if  ISSO  on  card  S-l 

is  not  0. 

Variable 

Format 

Column 

Description 

Z1 

E  10.  C 

1-10 

Z-coordinate  for  the  first  point  with  stress 
output. 

Y1 

E  10.6 

11  -20 

Y' -coordinate  for  the  same  point. 

(This  card  continued  next  page.) 
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Variable 

Format 

Column 

22 

E  10.6 

21-30 

Y2 

E  10.6 

31-40 

Z3 

E  10.6 

41-50 

Y3 

E  10.6 

51-60 

Z4 

E  10.6 

61-70 

Y4 

E  10.6 

71-80 

Note: 


Description 


S-  me  information  for  other  points  with 
stress  output.  Read  as  many  pairs  of 
coordinates  as  indicated  by  ISSO  on  the 
S-l  card. 


Repeat  above  data  for  each  stringer  type,  J  =  1,  ITSN. 


For  nonlinear  analysis,  go  to  the  P- IB  card;  for  linear  stress  analysis,  go  to  the 
0-1  card. 
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ANALYSIS  CONTROL 


P-1A  Parameter  Card:  Used  for  bifurcation  analysis,  INDIC  =  1  on  G-l  card. 


Variable 

Format 

Column? 

Description 

DELBIF 

E  10.6 

1-10 

Error  tolerance  in  power -Iteration  for 
eigenvalue.  If  DELBIF  is  zero  or  blank, 
the  error  tolerance  .  0001  is  used. 

SHIFT 

E  10.6 

11-20 

Initial  eigenvalue  shift,  if  any. 

IBOND 

15 

21-25 

1  -  The  boundary  conditions  for  incre¬ 
mental  displacements  are  different 
from  the  prebuckling  displacements. 

0  -  The  boundary  conditions  are  the  same 
for  incremental  and  prebuckling 
displacements. 

ISHIFT 

15 

26-30 

Number  of  eigenvalue  shifts  permitted. 

ITERAT 

15 

31-35 

Number  of  inverse  power  iterations 
permitted  between  shifts. 

If  IBOND  =  0,  go  the  0-1  card. 


P-1A1  Parameter  Card:  Incremental  displacement  boundary  condition.  Used  for  bi¬ 
furcation  analysis  with  different  boundary  conditions  for  incremental  and  prebuckling 
displacement  (INDIC  =  1  on  G-l  card  and  IBOND  =  l  on  P-1A  card). 


Boundary  lines  are  numbered  1  through  4  according  to  the  following  convention: 


Row  1 

-  Boundary  line  1 

Column  NC 

-  Boundary  line  2 

Row  NR 

-  Boundary  line  3 

Column  1 

-  Boundary  line  4 

Y  or  0,tj 


(This  card  continued  next  page . ) 
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Variable 


Format 


Columns  Description 

JBLN(^  415  1-20  Boundary  code  for  Line  I. 

0  —  Specified  on  P-1A2  cards 

1  -  Simple  support 

2  -  Clamped 

3  -  Unrestrained 

4  -  Symmetry 

5  -  Anti-symmetry 

6  -  Closed  shell 


Note:  See  discussion  of  B-l  card.  JBLN(I)  =  6  only  if  IBLN(I)  =  6. 


P2^2_^^amet^_C^^s:  Used  only  if  JBLN(I)  =  0  on  the  P-1A1  card.  One  P-1A2 
card  for  each  I  with  JBLN  =  0  .  The  cards  refer  to  boundary  lines  with  increasing 
numbers. 


Variable  Format  Columns 

JCOND(I)  4 15  1-20 

1=1,  4 


Description 

Freedom  of  movement  of  Line  I  in  regard 
to  the  w,  v,  u  ,  and  /?  displacements, 
where  w  ,  v  ,  and  u  are  displacements  in 
the  direction  of  the  coordinates  Z,  Y,  and 
X,  respectively  and  0  is  the  rotation 
around  a  tangent  to  the  edge. 

1  —  Free  to  move 
0  —  No  movement 


Y.n.v.e 


X,£,U 


Go  to  the  0-1  card. 
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| 

1 
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P-1D  Parameter  Card;  Used  for  nonlirear  analysis,  INDIC  =  2  or  3  on  G-l  card. 


Variable  Format  Columns  Description 

DELX  E  10.6  1-10  Error  tolerance. 

WUND  E10.6  11—20  Underrelaxation  factor. 

Note:  If  DELX  is  zero  or  blank,  the  error  tolerance  .  0001  Is  used.  If 

WUND  is  zero  or  blank,  the  relaxation  factor  of  1.0  is  used  initially. 
The  relaxation  factor  is  increased  internally  if  convergence  is  mono¬ 
tonic  but  slow,  and  it  is  decreased  internally  if  convergence  is  highly 
oscillatory.  If  WUND  *  0,  the  input  value  remains  unchanged 
regardless  of  convergence. 


Variable  Format  Columns  Description 

ISTART  15  21-25  Starting  Code:  (see  also  Section  7. 1) 

0  -  Begin  new  case. 

1  -  Restart  case  from  1st  record. 

2  -  Restart  case  from  2nd  record. 

3  -  Restart  from  3rd  record  Last  load  step. 

Note:  For  plasticity  (INDICT)  ISTART  is 
either  0  or  3.  ISTART  1  or  2  cannot 
be  used. 


ISEC 

15 

26-30 

Number  of  CPU  seconds  of  run  time  at 
which  run  should  be  terminated  and  data 
saved  on  restart  tape.  (See  Strategy, 
Section  4) 

ICUT 

15 

31-35 

Total  number  of  times  step  size  may 
be  cut.  (See  Strategy,  Section  4). 

INEWT 

15 

36-40 

Number  of  Newton  steps  which  may  be 
taken.  (See  Strategy,  Section  4) 

ISTRAT 

15 

41-45 

Number  of  times  step  size  is  cut  between 
Newton  steps.  (See  Strategy,  Section  4) 
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OUTPUT 


0-1  Output  Definition  Card 

Variable 

Format 

Columns 

Description 

IPR1 

15 

1-5 

Print  displacement  solutions  every  IPR1 
load  step.  * 

IPR2 

15 

6-10 

Print  stress  resultants  every  IPR2  load 
step.  * 

IPR3 

15 

11-15 

Frint  stresses  every  IPR3  load  step.* 

♦Note: 

With  IPR1,  IPR2,  or  IPR3  equal  zero  corresponding  output  is  suppressed. 
This  option  can  be  used  also  in  the  linear  or  bifurcation  buckling  case  if 
one  wants  to  suppress  all  regular 1  output  in  favor  of  the  selected  output 
defined  by  IPX,  IPY. 

Variable 

Format 

Columns 

Description 

IPX 

15 

16-20 

Number  of  selected  rows  along  which  dis¬ 
placements  and/or  stress  resultant 
solutions  will  be  printed.  IPX  £  50 

IPY 

15 

21-25 

Number  of  selected  rclumns  along  which 
displacements  and/or  stress  resultant 
solutions  will  be  printed.  IPY  <  50 

IPRD 

15 

26-30 

Print  selected  displacement  solutions 
every  IPRD  load  step. 

IPRS 

15 

31-35 

Print  selected  stress  resultant 
solutions  every  IPRS  load  step. 

IPLOT 

15 

36-40 

See  Section  8  for  details. 

0  -  Ho  plots. 

1  -  Plot  and  print  special  output  (see  notes 

below).  Plot  selected  output  (if  any). 

2  -  Plot  special  output.  Plot  selected 

output  (if  any). 

IHIST 

15 

41-45 

See  Section  8  for  details. 

The  number  of  selected  points  for  which 
history  plots  are  required.  If  no  history 
plots  are  required,  set  IHIST=0. 

IHIST  <4. 

(This  card  continued  next  page . ) 
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Notes:  1.  Plot  means  data  will  be  saved  on  special  tape  for  use  with 
postprocessor. 

2 .  Special  output; 

If  INDICT  (on  the  G-l  card)  -  All  displacements  and  stress  resultants. 
If  INDIOl  -  All  prebuckling  displacements  and  stress  resultants  and 
all  displacements  in  the  buckling  mode. 

If  INDIC-2  or  3  -  All  displacements  and  stress  resultants  for  the  last 
load-step  and  the  difference  between  the  last  two  dis¬ 
placement  solutions  (collapse  mode). 

3.  For  IPLOT=l  or  2  and  INDIC-2  or  3,  a  history  plot  for  the  maximum 
displacement  is  provided  in  addition  to  IHIST  plots. 

4.  An  effort  was  made  to  eliminate  duplication  of  output  printing. 


0-2  Selected  Output  Along  Rows 


Variable  Format  Columns  Description 

XPL(I)  8E  10.6  1-80  X  coordinates  of  rows  along  which 

I  =  1,  IPX  selected  output  is  desired. 

Note:  Include  these  cards  only  if  IPX  on  card  0-1  greater  than  zero  (selected 
output  is  requested).  Card  format  is  8E  10.6,  therefore  more  than  one 
0-2  card  might  be  necessary  to  input  all  X  coordinates. 


0-3  Selected  Output  Along  Columns 


Variable  Format  Columns  Description 

YPL(I)  8E  10.6  1-80  Y  coordinates  of  columns  along  which 

1  =  1,  IPY  selected  output  is  desired. 


Note:  Include  these  cards  only  if  IPY  on  card  0-1  greater  than  zero  (selected 
output  is  requested).  Card  format  is  8E  10.6,  therefore  more  than  one 
0-3  card  might  be  necessary  to  input  all  Y  coordinates. 
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0-4  Selected  Points  for  History  Plot 


ThiB  card  should  be  included  only  if  IHIST  on  card  0-1  is  greater  than  zero. 
The  number  of  entries  on  this  card  will  be  equal  to  2  x  IHIST. 


Variable 

Format 

Columns 

Description 

XI 

E10.6 

1-10 

X  -  Coor  llnatt  for  first  mes  l  point  at 
which  history  data  are  paved. 

Y1 

E10.6 

11-20 

Y  -  Coordinate  for  first  u.-sh  point  at 
which  history  data  are  saved. 

X2 

E10.6 

21-30 

Y2 

E10.6 

31-40 

X3 

E10.6 

41-50 

Y3 

E10.6 

51-60 

X4 

E10.6 

61-70 

Y4 

E10.6 

71-80 

If  the  plasticity  branch  is  used  there  are  no  additional  cards  to  be  read. 
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WALL  CONSTRUCTION 


M-l  Wall  Type  Card 


Cards  on  which  data  are  to  be  read  for  the  different  types  of  shell  wall  are  given  in 
parentheses. 


Variable  Format 
IWALL  15 


Columns  Description 

1-5  1-A  user-written  subroutine,  WALL,  is 

going  to  be  provided.  (Except  for  the 
options  under  IWALL  -  8  or  9,  this  is 
the  only  way  to  introduce  variable 
properties).  Notice  that  the  subroutines 
corresponding  to  IWALL  =  2  through  7 
may  be  cal’.  1  from  WALL.  Input  is 
complete. 

2—  Monocoque  shell  (may  be  orthotropic). 
(M-2B) 

3—  Skew-stiffened  shell.  (M-2C1,  M-2C2) 

4—  Layered,  fiber-wound  shell.  (M-2D1, 
M-2D2) 

5—  Layered  shells,  layers  may  be 
orthotropic.  (M-2E1,  M-2E2) 

6  — Corrugated  shell.  (M-2F) 

7—  Corrugated  shell  with  smooth  skin. 
(M-2G1,  M-2G2) 

8—  Orthotropic  shell.  (Temperature  and 
wall  properties  are  defined  by  user 
written  MATER  subroutine  and  may 
vary  through  the  thickness  as  well  as 
with  shell  coordinates.)  No  more  cards 
to  read  unless  NSTRI  or  NRING  below 
equals  1. 

9—  Orthotropic  shell.  (Temperature  and 
wall  properties  are  read  at  selected 
mesh  points  and  may  vary  through 
the  thickness  as  well  as  with  shell 
coordinates. )(V-1,  V-2,  V-3,  V-4) 


(This  card  continued  next  page . ) 
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Variable 

Format 

Columns 

Descriotion 

NSTRI 

15 

6-10 

0-  No  smeared  stringers. 

1  -  Smeared  stringers  (input  data  on 

card  N-l). 

NRING 

15 

11-15 

0  -  No  smeared  rings 

1  -  Smeared  rings  (input  data  on  card 
N-3). 

Unless  IWALL  is  8  or  9,  this  card  is  complete. 

Discrete  stiffeners  can  be  added  as  desired  above  (S-l  through  S-7). 


Variable 

Format 

Columns 

Description 

IP 

15 

16-20 

Number  of  points  across  wall  thickness. 
Must  be  odd  number  and  not  less  than  3  or 
more  than  9.  Omit  unless  IWALL  =  8  or  9. 

IM 

15 

21-25 

Number  of  rows  selected  for  input  of 
temperature  or  properties.  Omit  unless 
IWALL  =  9.  IM  >  2. 

JM 

15 

26-30 

Number  of  columns  selected  for  input  of 
temperature  or  properties.  Omit  unless 
IWALL  =  9.  JM>  2. 
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MATERIAL  PROPERTIES 


M-2B  (IWALL  =  2,  Monocoque  Shell) 


Note:  If  shell  is  isotropic,  punch  only  the  first  four  fields  on  this  card. 
For  an  orthotropic  shell  with  EY1  =  0,  punch  a  small  but  nonzero 


value  for  EY1. 


Variable 

Format 

Columns 

AT 

E  10.6 

1-10 

EX1 

E  10.6 

11-20 

XNU 

E  10.6 

21-30 

Z 

E  10.6 

31-40 

EY1 

E  10.6 

41-50 

G 

E  10.6 

51-60 

M-2C1  (IWALL 

=  3,  skew-stiffened  shell) 

Variable 

Format 

Columns 

T3 

E  10.6 

1-10 

E3 

E  10.6 

11-20 

U3 

E  10.6 

21-30 

TH 

E  10.6 

31-40 

A 

E  10.  G 

41-50 

B 

E  10.6 

51-60 

H3 

E  10.6 

61-70 

AK3 

E  10.6 

71-80 

Description 
Wall  thickness. 

Young's  Modulus  in  X-direction. 

Poisson's  Patio  -  nxy  ,  (Mxy  Ex  =  #iyx  Ey) 

Distance  from  reference  surface  to 
shell  midsurface  (positive  if  the  mid¬ 
surface  is  outside  of  the  reference 
surface,  i.e.,  /.coordinate  for  mid¬ 
surface  is  positive). 

Young's  Modulus  in  Y-direction  (however 
shell  is  assumed  Isotropic  if  EY1  -  0). 

Shear  Modulus. 


Description 

Wall  thickness 
Young's  modulus 
Poisson's  ratio 

Angle  (deg)  between  stiffeners  and 
X-coordinatc  lines 

Stiffener  spacing  (along  Y-cooidinatc 


M-2C2  (RVALL  =  3,  skew- stiffened  shell) 


Variable 

Format 

Columns 

Description 

Z 

E  10.6 

1-10 

Distance  from  reference  surface  to  skin 
midsurface  (positive  if  the  midsurface 
is  outside  the  reference  surface,  i.e., 
at  positive  Z- coordinate). 

M-2D1  (WALL  =  4,  layered,  fiber- wound  shell) 

Variable 

Format 

Columns 

Description 

EF 

E  10.6 

1-10 

Young's  modulus  of  fibers 

EM 

E  10.6 

11-20 

Young's  modulus  of  matrix 

UF 

E  10.6 

21-30 

Poisson's  ratio  of  fibers 

UM 

E  10.6 

31-40 

Poisson's  ratio  of  matrix 

Z 

E  10.6 

41-50 

Distance  from  reference  surface  to  shell 
midsurface  (positive  if  midsurface 
is  outside  the  reference  surface;  i.e., 
at  a  positive  Z  coordinate) 

LAYERS 

15 

51-55 

Number  of  layers.  LAYERS  <  20. 

M-2D2  (WALL  =  4,  layered,  fiber-wound  shell) 

Variable 

Format 

Columns 

Description 

TT(J) 

E  10.6 

1-10 

Thickness  of  layer  (inner  layer  has 
index  1 ,  outer  layer  has  index  LAYERS) 

XX  (J) 

E  10.6 

11-20 

Matrix  content  (by  volume)  in  percent/100 

BE(J) 

E  10.6 

21-30 

Winding  angle 

0(J) 

E  10.6 

31-40 

Contiguity  factor 

Note: 

Repeat  Card  M-2D2  for  J  =  1, 

LAYERS. 

M-2E1  (IWALL  =  5, 

layered  orthotropic  shells) 

Variable 

Format 

Columns 

Descripi'on 

Z 

E  10.6 

1-10 

Distance  from  reference  surface  to  shell 
midsurface  (positive  if  midsurface  is 
outside  of  reference  surface,  i.e. ,  at  a 
positive  Z  coordinate). 

LAYS 

15 

11-15 

Number  of  layers.  LAYS  <  20. 
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M-2K2  (IWALL  =  5,  layered  orthotropic  shell) 


Variable 

Format 

Columns 

Description 

TL(J) 

E  10.6 

1-10 

Layer  thickness  (inner  layer  No.  1,  outer 
layer  No.  LAYS) 

EX5(J) 

E  10.6 

11-20 

Modulus  of  elasticity  in  x  direction 

EY5(J) 

E  10.6 

21-30 

Modulus  of  elasticity  in  y  direction 

UXY(J) 

E  10.6 

31-40 

Poisson's  ratio  -  u  ,  (u  E  =u  E  1 
^xy’  ,fxy  x  Myx  y' 

G5(J) 

E  10.6 

41-50 

Shear  modulus 

Note:  Repeat  Card  M-2E2  for  J  = 

1.  LAYS. 

M-2F  (IWALL 

=  6,  Corrugated  Shell) 

Variable 

Format 

Columns 

Description 

CT6 

E  10.6 

1-10 

Thickness  of  corrugated  sheet 

E6 

E  10.6 

11-20 

Young's  modulus 

U6 

E  10.6 

21-30 

Poiss9n's  ratio 

CC6 

E  10.6 

31-40 

CH6 

E  10.6 

41-50 

See  figure 

CD6 

E  10.6 

51-60 

CB6  E  10.6  61-70  Centerllne-to-centerline  spacing  of 

corrugations 


(This  card  continued  next  page . ) 
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Variable  Format  Columns 


Description 


Z  E  10. 6  71-80  Distance  from  reference  surface  to 

shell  midsurface  (positive  if  midsurface 
of  shell  is  outside  the  reference  surface, 
i.e. ,  at  a  positive  Z  coordinate) 


M-2G1  (IWALL  =  7,  Corrugated  Shell  with  one  smooth  skin) 


Variable 

Format 

Columns 

Description 

CT7 

E  10.6 

1-10 

Thickness  of  corrugated  sheet 

E7 

E  10.6 

11-20 

Young's  modulus 

U7 

E  10.6 

21-30 

Poisson's  ratio 

CC7 

E  10.6 

31-40 

CH7 

E  10.6 

41-50 

See  figure 

CD7 

E  10.6 

51-60 

4 

CB7 

E  10.6 

61-70 

C  enterline- to-centerline  spacing  of 
corrugations 

Z 

E  10.6 

71-80 

Distance  from  reference  surface  to  mid¬ 
surface  of  shell  (positive  if  midsurface 
is  outside  of  reference  surface;  at 
positive  Z  coordinate) . 
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M-2G2  (IWALL  =  7,  Corrugated  Shell  with  one  smooth  skin) 


Variable 

Format 

Columns 

Description 

TS 

E  10.6 

1-10 

Thickness  of  smooth  skin 

ES 

E  10.6 

11-20 

Young's  modulus  of  skin 

US 

E  10.6 

21-30 

Poisson's  ratio  of  skin 

PHI 

E  10.6 

31-40 

Reduction  factor  for  torsional  stiffness 

ANC 

E  10.6 

41-50 

0.-  Inside  corrugation 

1.-  Outside  corrugation 

SMEARED  STIFFENERS 


Note:  For  smeared  rectangular  stiffeners,  GJ is  computed  by  the  program  as 
H  x  T3/3  when  H  >  T  (somewhat  unconservative  for  H  T). 

N-l  Smeared  Stringers 


Read  only  if  NSTRI  =  1  on  Card  M-l. 


Variable 

Format 

Columns 

El 

E  10.6 

1-10 

U1 

E  10.6 

11-20 

011 

E  10.6 

21-30 

D1 

E  10.6 

31-40 

AK1 

E  10.6 

41-50 

N-2A  Smeared  Rectangular  Stringers 

Read  only 

if  011  =  0  on  Card 

N-l. 

Variable 

Format 

Columns 

T1 

E  10.6 

1-10 

HI 

E  10.6 

11-20 

N-2B  Smeared  Arbitrary  Stringers 

Read  only 

if  011  =  1  on  Card  N-l. 

Variable 

Format 

Columns 

A1 

E  10.6 

1-10 

SI1 

E  10.6 

11-20 

XII 

E  10.6 

21-30 

SI 

E  10.6 

31-40 

EZ1 

E  10.6 

41-50 

HI 

E  10.6 

51-60 

Description 

Young's  modulus 

Poisson's  ratio 

Ol  -  Rectangular  stringers 
1.-  Arbitrary  stringers 

Stringer  spacing  (arc  length) 

0.-  Internal  stringers 
1.-  External  Btringers 


Desc  rjption 

Stringer  thickness 
Stringer  height 


Description 
Area  of  stringer 

Moment  of  inertia  of  stringer  about  V-axis 

Moment  of  inertia  of  stringer  about  Z'-axis 

Torsional  stiffness  of  stringer,  GJ 

Eccentricity  of  stringer.  (See  figure 
for  discrete  stringers.) 

Stringer  Height 
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N-3  Smeared  Rings 


Head  only  If  NRING  =  1  on 

Card  M-i. 

Variable 

Format 

Columns 

Description 

E2 

E  10.6 

1-10 

Young's  Modulus 

U2 

E  10.6 

11-20 

Poisson's  Ratio 

012 

E  10.6 

21-30 

0.-  Rectangular  rings 

1.—  Arbitrary  rings 

D2 

E  10.6 

31-40 

Ring  spacing  (arc  length) 

AK2 

E  10.6 

41-50 

0.-  Internal  rings 

1.—  External  rings 

N-4A  Smeared  Rectangular  Rings 

Read  only  if  012 

=  0  on  Card  N-3. 

Variable 

Format 

Columns 

Description 

T2 

E  10.6 

1-10 

Ring  Thickness 

H2 

E  10.6 

11-20 

Ring  Height 

N-4B  Smeared  Arbitrary  Rings 

Read  only  if  012 

=  1  on  Card  N-3. 

Variable 

Format 

Columns 

Description 

A2 

E  10.6 

1-10 

Area  of  ring 

SI2 

E  10.6 

11-20 

Moment  of  inertia  of  ring  about  X'-axis 

XI2 

E  10.6 

21-30 

Moment  of  inertia  of  ring  about  Z'-axis 

S2 

E  10.6 

31-40 

Torsional  stiffness  of  ring,  GJ 

EZ2 

E  10.6 

41-50 

Eccentricity  of  ring  (See  figure  for  dis¬ 
crete  rings) 

H2 

E  10.6 

51-60 

Ring  Height 

Note:  Stringers  run  parallel  with  X  and  rings  parallel  with  Y  coordinate  axis. 

The  remaining  cards  are  to  be  read  only  if  IWALL  =  9. 
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VARIABLE  WALL  PROPERTIES 


V-I  (For  Orthotropic  Shell  with  variable  properties  -  IWALL  =  9) 


Variable 

Format 

Columns 

Description 

TD 

E  10.6 

1-10 

Wall  thickness 

Z 

E  10.6 

11-20 

Distance  from  reference  surface  to  mid- 
suriace  of  shell  (positive  if  midsurface  is 
outside  of  reference  surface,  i.e. ,  at  a 
positive  Z  coordinate) 

V-2  X-Coarse  Grid  Cards 

Variable 

Format 

Columns 

Description 

XM(I) 

I  =  1,IM 

8E  10.6 

1-80 

X-coordinates  of  selected  mesh  (for  IM, 
see  M-l  card) 

Note: 

Card  format  is  8E  10.6;  therefore,  more  than  one  V-2  card  might  be 
necessary  to  input  all  X  coordinates. 

V-3  Y-Coarse  Grid  Cards 

Variable 

Format 

Columns 

Description 

YM(J) 

J  =  1,  JM 

8E  10.6 

1-80 

Y-coordinates  of  selected  mesh  (for  JM, 
see  M-l  card) 

Note:  Card  format  is  8E  10.6;  therefore, more  than  one  V-3  card  might  be 

necessary  to  input  all  Y  coordinates. 


The  following  cards  are  reed  in  a  loop  over  the  selected  mesh  points.  They  are  read 
either  by  row  or  by  column,  depending  on  which  direction  has  the  largest  number  of 
original  gridpoints.  With  NR  (the  number  of  points  in  the  X-direction),  NC  (the  number 
of  points  in  the  Y -direction  -  see  Card  D-I),  and  IM  and  JM  (the  number  of  selected 
points  in  these  directions) 
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for  NC  i=  NR 

for  NC  >  NR 


MSTA  JM 
NSTA  IM 

MSTA  =  IM 
NSTA  =  JM 


The  data  card  is  repeated  for  L  =  1,  IP  (number  of  points,  numbered  inner  to  outer, 
see  M-l  card) 


M  =  1,  MSTA 
N  =  1,  NSTA 

Note  that  L  is  the  inner  loop,  M  is  the  middle  loop,  and  N  is  the  outer  loop. 


V-4  Variable  Properties  Cards 


Variable 

Format 

Columns 

Description 

TDEG(L,  M,  N) 

E  10.6 

1-10 

Wall  temperature 

EX(L,  M,  N) 

E  10.6 

11-20 

Modijlus  of  elasticity  in  X  direction 

EY(L,  M,  N) 

E  10.6 

21-30 

Modulus  of  elasticity  in  Y  direction 

U(L,M,N) 

E  10.6 

31-40 

Poisson’s  ratio  -  u  .  Ui  E  -a  E 
xy’  ''xy  x  yx  ; 

G(L,  M,N) 

E  10.6 

41-50 

Shear  modulus 

Al(L.M.N) 

E  10.6 

51-60 

Coefficient  of  thermal  expansion  in 

X  direction 

A2(L,  M,  N) 

E  10.6 

61-70 

Coefficient  of  thermal  expansion  in 

Y  direction 

Note:  Card  format  is  7E  10.  6,  therefore  many  V-4  cards  will  be  necessary 
to  input  all  wall  properties. 


END  OF  INPUT 
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Section  7 

USE  OF  STAGS  PROGRAM 


This  section  illustrates  the  use  of  the  STAGS  computer  program  to  solve  some  typical 
examples.  It  also  illustrates  the  form  and  the  interpretation  of  STAGS-generated  out¬ 
put.  To  conserve  space,  only  representative  portions  of  the  output  are  shown.  Cases 
have  been  chosen  which  are  deemed  suitable  for  program  checkout.  Hence ,  together 
they  cover  most  of  the  branches  of  the  program ,  and  they  are  relatively  inexpensive  to 
run.  It  is  not  intended  to  demonstrate  the  capability  of  the  program;  many  of  the  cases 
can  be  solved  by  simpler  means. 

7. 1  SAMPLE  CASE  1  -  CYLINDRICAL  SHELL  SEGMENT 

Consider  a  circular  cylindrical  shell  panel  subjected  to  a  radial  load  at  the  center 
(Fig.  7-1).  The  panel  is  simply  supported  along  its  curved  boundaries  and  free  along 
its  straight  boundaries.  Because  of  the  symmetrical  arrangement,  only  one  quarter 
of  the  panel  need  be  analyzed.  This  is  indicated  in  the  figure;  the  boundaries  are  iden¬ 
tified  as  shown.  Figure  7-2  shows  the  input  cards  associated  with  this  case;  Fig.  7-3 
shows  portions  of  the  output. 

Computer  output  begins  with  a  page  containing  pertinent  input  data  and  stiffness  co¬ 
efficients  followed  by  the  determinant  of  the  system  of  equations  and  an  indication  of 
the  elapsed  time.  The  linear  solution  is  printed  in  five  segments  (according  to  the 
0-1  type  input  card):  first,  the  displacement  field;  second,  the  moment  and  stress  re¬ 
sultants  field;  third,  the  stress  field;  fourth,  selected  displacements;  fifth,  selected 
moment  and  stress  resultants.  For  cases  in  which  the  shell  is  loaded  by  a  fixed  dis¬ 
placement  (unit  end  shortening),  it  is  of  interest  to  find  the  applied  load.  Therefore, 
the  stress  resultants  are  integrated  at  the  NR-1  row  and  the  result  is  given  as  output. 
In  this  particular  case,  the  axial  load,  near  zero,  is  of  no  interest. 
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After  the  linear  solution,  the  program  begins  to  compute  the  nonlinear  displacements 
for  successively  increasing  load  factors  as  specified  by  means  of  the  load  increment 
parameter.  Only  the  selected  output  related  to  load  steps  1,  2,  and  6  is  shown.  For 
each  load  step,  the  load  factors  PA  and  PB  for  load  system  A  and  load  system  B  are 
printed.  The  load  factor  PA  for  this  last  step  was  6.0 ,  and  the  solution  required  1 
iteration.  For  each  iteration ,  the  following  information  is  displayed:  maximum  dis¬ 
placement  change,  displacement  component,  location,  relative  error,  and  relaxation 
factor. 

The  nonlinear  solution  is  displayed  either  in  its  entirety  as  in  the  linear  solution  or 
only  at  selected  mesh  points  as  requested  by  the  user  on  the  0-1  type  of  input  card. 

When  program  termination  occurs  because  of  an  internal  input  card  parameter  (see 
L-l  and  P-1B  input  cards),  a  statement  describing  the  reason  for  termination  will 
appear  in  the  output  and  the  last  three  successive  solutions  will  be  saved  on  tape  18 
(output  tape).  A  second  run  may  start  from  any  of  these  three  records,  provided 
ohanges  are  made  to  the  input  data,  such  as  starting  load  factor  (L-l  type  input  card) 
and  restart  command  (P-1B  type  input  card). 

A  display  of  input  cards  and  selected  output  for  the  second  run  follows  the  first  run 
(Figs.  7- 3a  and  7 -3b). 

Note  that  the  load  step  increment  may  be  cut  or  increased  automatically  based  on  an 
internal  criteria. 

The  analysis  may  be  continued  until  the  collapse  load  has  been  determined ,  or  con¬ 
vergence  difficulties  for  small  load  increments  are  encountered. 
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Fig.  7-1  Sample  Case  1  -  Cylindrical  Shell  Panel 
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Fig.  7-2  Display  of  Input  Cards  for  Sample  Case  1,  First  Run 
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Fig.  7-3  Excerpt  of  Output  for  Sample  Case  1,  First  Run. 
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Fig.  7-3  (Cont.) 
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The  second  nonlinear  run  is  a  continuation  of  the  first  run,  which  terminated  (because 
the  maximum  load  was  obtained)  at  load  step  6  with  PA  equal  G.  0.  Because  the  run 
was  started  with  the  3rd  record,  i.  e.  ISTART  equal  3  (see  p-lB  type  input  card),  the 
starting  load  factor  must  be  changed  accordingly.  Hence,  STLD  is  set  equal  6. 0  (see 
L-l  type  input  card). 


SAMPLE  CASE  1 
SECOND  RUN  -  INPUT 
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Fig.  7~3a  Display  of  Input  Cards  for  Sample  Case  1,  Second  Run 
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SAMPLE  CASE  1,  SECOND  RUN  -  OUTPUT 
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Fig.  7-3b  Excerpt  of  Output  for  Sample  Case  1,  Second  Run 
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Fig.  7-3b  (Cont.) 


7.2  SAMPLE  CASE  2  -  CYLINDRICAL  SHELL  WITH  ELLIPTIC  CROSS  SECTION 


This  case  demonstrates  calculation  of  the  buckling  load,  according  to  bifurcation 
theory,  for  a  cylinder  with  elliptic  cross  section  (Fig.  7-4)  and  subjected  to  the  combined 
effects  of  a  line  load  along  boundary  line  1  (load  pattern  A)  and  uniform  internal  pressure 
(load  pattern  B).  Here  will  be  determined  the  critical  axial  load  corresponding  to  an 
internal  pressure  of  20  psi.  Because  of  symmetry  ,  only  one  eighth  of  the  cylinder 
needs  to  be  analyzed.  In  the  prebuckling  analysis,  the  restricting  assumption  is  made 
that  the  buckling  pattern  also  is  symmetric  about  the  same  two  planes.  The  input  cards 
associated  with  this  case  are  illustrated  in  Fig.  7-5.  Portions  of  the  output  are  dis¬ 
played  in  Fig.  7-6. 

After  displaying  the  input  parameters,  stiffness  coefficients,  and  the  determinant  of 
the  system,  the  linear  solution  is  given  as  in  the  previous  example,  followed  by  the 
determinant  of  the  system  based  on  the  new  boundary  conditions  (if  the  incremental  and 
prebuckling  displacement  boundary  conditions  are  different)  and  an  initial  shift  if  any. 

The  eigenvalue  (buckling  load)  is  computed  by  a  series  of  iterations  displaying:  itera¬ 
tion  count,  Rayleigh  quotient  and  the  accelerated  estimate  of  the  eigenvalue.  To  re¬ 
duce  the  number  of  iterations  required  for  convergence  an  initial  shift  (lower  than  the 
expected  solution)  was  utilized. 

As  a  final  result,  the  normalized  buckling  mode  is  printed  in  the  same  format  as  the 
linear  solution. 
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Fig.  7-4  Sample  Case  2  -  Elliptic  Cylinder 


SAMPLE  CASE  2 
INPUT 
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Fig.  7-5  Display  of  Input  Cards  for  Sample  Case  2 
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SAMPLE  CASE  2  -  OUTPUT 
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7-«  Excerpt  of  Output  for  Sample  Ca*e  2 
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7.3  SAMPLE  CASE  3  -  CYLINDER  WITH  RECTANGULAR  CUTOUT 


This  case  demonstrates  calculation  of  the  buckling  load,  according  to  bifurcation  theory, 
for  a  cylinder  with  rectangular  cutout  (Fig.  7-7)  and  subjected  to  a  uniform  line  load 
applied  on  boundary  line  1.  Because  of  symmetry,  only  one  eighth  of  the  shell  needs 
to  be  analyzed.  The  input  cards  associated  with  this  case  are  illustrated  in  Fig.  7-8. 
Portions  of  the  output  are  presented  in  Fig.  7-9.  The  output  is  the  same  as  in  Sample 
Case  2. 
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LINE  3 

Fig.  7-7  Sample  Case  3  -  Cylinder  With  Rectangular  Cutout 

SAMPLE  CASE  3 
INPUT 
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SAMPLE  CASE  3 
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Fig.  7-8  Display  of  Input  Cards  for  Sample  Case  3 


7-23 


5 


Fig.  7-9  Excerpt  of  Output  for  Sample  Case  3 
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Fig.  7-9  (Cont.) 
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Fig.  7-9  (Cont. 
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Fig.  7-9  (Cent.) 


7.4  SAMPLE  CASE  4  -  TOROIDAL  SHELL 


The  critical  axial  load  according  to  bifurcation  theory  is  to  be  determined  for  a  toroidal 
shell  segment  simply  supported  at  circular  edges  and  stiffened  with  internal  stringers. 
Because  this  is  a  shell  of  revolution  with  axially  symmetrical  loading,  it  could  be 
analyzed  with  a  one  dimensional  computer  program.  Such  cases  are  included  here 
because  they  provide  an  opportunity  to  check  the  program.  The  geometry  of  the  shell 
is  shown  in  Fig.  7-10. 

It  is  assumed  here  that  the  critical  axial  load  is  obtained  for  a  mode  with  36  circumfer¬ 
ential  waves.  Therefore,  a  shell  covering  one  half  the  length  and  1/144  of  the  circum¬ 
ference  can  be  analyzed.  That  is,  the  shell  segment  is  expected  to  buckle  with  one 
quarter  of  a  sine  wave  in  the  circumferential  direction.  Thus,  on  the  sides  formed  by 
the  generators  (2  and  4),  symmetry  will  be  used  for  both  in  the  prebuckling  analysis 
and  symmetry  for  one  and  antisymmetry  for  the  other  in  the  buckling  analysis.  The 
input  cards  associated  with  this  case  are  displayed  in  Fig.  7-11.  Portions  of  the  output 
are  presented  in  Fig.  7-12. 

For  zero  initial  eigenvalue  shift,  the  buckling  load  obtained  was  -1780  Ib/in.  (i.e. 
tension).  To  obtain  a  compressive  buckling  load,  a  run  was  performed  with  an  initial 
eigenvalue  shift  of  +1780,  but  the  results  were  again  -1780  lb/in.  This  indicates  that 
the  compressive  buckling  load  must  be  greater  than  5340.  Therefore,  an  initial  eigen¬ 
value  shift  of  5340  was  tried  leading  to  a  bifurcation  load  of  6390  lb/in. 
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Fig.  7-10  Sample  Case  4  -  Toroidal  Shell 
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Fig.  7-11  Display  of  Input  Cards  for  Sample  Case  4 
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7.5  SAMPLE  CASE  5  -  CORRUGATED  CYLINDER 


The  critical  axial  line  load  according  to  bifurcation  theory  is  to  be  determined  for  a 
simply  supported  corrugated  circular  cylindrical  shell  stiffened  with  internal  rings. 
Because  this  is  a  shell  of  revolution  with  axially  symmetrical  loading,  it  could  be 
analyzed  with  the  BOSOR4  computer  program  (Ref.  15)  leading  to  a  critical  axial 
line  load  of  952  lb/in.  for  a  mode  with  7  circumferential  waves.  Therefore,  a  shell 
segment  covering  one  half  the  length  and  1/28  of  the  circumference  can  be  analyzed. 

That  is,  the  shell  segment  is  expected  to  buckle  with  one  quarter  of  a  sine  wave  in  the 
circumferential  direction.  Thus,  the  boundary  condition  on  lines  2  and  4  is  identical 
to  that  of  Sample  Case  4  both  for  the  prebuckling  and  for  the  buckling  analysis.  The 
geometry  of  the  shell  is  shown  in  Fig.  7-13.  The  input  cards  associated  with  this 
case  are  displayed  in  Fig.  7-14.  Portions  of  the  output  are  presented  in  Fig.  7-15. 

To  check  the  branch  for  arbitrary  stiffeners,  the  rectangular  rings  with  4.0-in.  spacing 
were  described  in  the  input  as  arbitrary  rings,  and  to  reduce  the  number  of  iterations, 
an  initial  eigenvalue  shift  at  1000.0  was  utilized.  The  critical  axial  line  load  determined 
here  is  987  lb/in. 
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SAMPLE  CASE  5  -  OUTPUT 


O 

u 

u 


o 

o 

o 


0) 

CO 

cd 

U 

0) 


8 

o 

fj 


o 

*  • 
s 


£  o 

WS 

y 


cd 

CO 

o 


u 

8 


o 

u 


•f  o  o  o 

ri  ar  *  «r 


JCXX 
ft  XII 

t  ►  ►> 

Mini/iw 


y  uj  y  Uf 


*  •  o 

M 

CftOUK 

Wt  MM  k 

U  J  V  K  tu 
m  ¥-  *-  i: 
JW  W  W  I 

l  J  mm 

nai.r*- 
«t  ►  ►  z 

VtVI  « 
13 

rirtv'nw 

HMHHH 

SImn«* 

Syyyy 

mmmh 
ft  -i  -»  -i  -J 


D 

or 

ft 


«  u 
J  u 
90 
O 


w  o 

8° 


19 

O 

£ 

UJ 

t 

►4 


* 

X 

o 


fr 

0 

o 

K 

W 


I 

t> 

U) 

E 


•  4  4  <  4 


►*  X  »- 


VI  w  M 


2  w  E 


W  »  mm* 


o 

u 


fe 

9 


£?92S 

onuo 

*  o  o  u  u 

r* 

I  ►  ►  ►  ► 

ft  ft  u-  a 

4  4  *1  4 

■  c.uao 
m  /  r  z 

*  :»  :»  ^ 
ft*  r»  f>  i  *  n 

/cr  ruo 


2522? 

r  o  o  c»  u 

o  u  o  o  o 

►  ►  ►  ►  ► 
ft  ft  ft  ft  ft 

4  4  4  4  4 

o  onii  ti 

5  S  5  5  § 


d 

U. 

© 

* 

D 

-> 

8 


K  II 
l>  M 


O  ~ 
-J  O 
J  u 
O  <1 


4  m  J 


r«c  rnnottiw  stiffness  coefficients  »«  CM.cot.mo  in  sushoutine  stiff 

CCCd.ll  CCCII.EI  CCC(I.I)  CCCtl.N)  cecti.si  cccti.tt 


58 

to  X 


c- 

Ui 

X 


O  Ui 
X  x 

to  to 


D  •- 
K.  VI 

X  UI 
O  3 
U.  O 
toi 

toi  ar 

o 

x  o 

toi  m 

to 

W  to. 
to.  O 


*  toi 
o  n 


O  H 

O  to* 
Ui  tot 

te  » 


o  o 

to.  to 


V* 


O  Ul 
'J 

n  o 

«*  toi 

•c  v 
«  o 

X  M 


V  to « 


2w 
j  »* 


h 

■ 

7  VI 
O  o 


<1  o 
s.  w 
to  vi 


•j  ■* 

m  • 

*-  to 


o  « 
v* 
►  C 

-J  -T 

e 

T  M 
Ui  W 

*r  vs 


o  to 
•  3 

ui  r 


M  «  O  O 
N  D  UJ  m 

to  o  *- 

^  to<  Ui  to. 
«  J  O 
x  -i  ft 
to  to 

v>  «  o  z 


T  • 

O  «J 

*-♦  to 


v*  o 

to  ui  to 

°SS 

1-70" 
7  Ui  to 

nr 
2  x  2 
•  «N  tot* 
X  MO 


• 

• 

• 

• 

• 

to 

«4 

to 

to 

to 

n 

K 

to 

K 

to 

to 

o 

4- 

O 

to 

w* 

o 

*4 

to 

to 

UI 

o 

fN 

«4 

C 

•4 

Ui 

• 

■ 

> 

■ 

to 

K  to  to 

W 

Ui 

N 

*n 

to 

M  to  FN 

Ui 

UI 

Ui 

to 

Ui 

1 

1  1  • 

•  X 

a 

a 

O 

a 

~  Ui 

Ui  Ui  Ui 

O  4 

Ui 

to 

• 

O  IN 

m  to  4"* 

4- 

4- 

X 

4- 

0 

2  IN 

o  to  ^ 

»» 

M  K 

to  to  ^ 

So 

o 

o 

8 

O 

sr 

•  r  r 

4  •  •  • 

a  o 

Ui 

Ui 

o 

Ui 

to  •«  o  r  -4  f>y  ( 

X  to 

to 

to 

te  to 

to 

_ 

X  1 

1  • 

O  3 

3 

3 

Ui 

3 

■ 

o 

m 

44 

to 

to 

to 

X 

to 

X 

to 

to  o 

o 

o 

3 

o 

Ui 

Utototototototo 

rttDCBOVB 

Dtototototototo 


3  uieeooeo 
ii«o  ♦  ♦  •  ♦  ♦  ♦ 
totowutuwwut 

t  Ut  W.  to  X)  to  to  to  to 

r  or  of  c  ®  «■  B  k  c 
m  D  to  to  v  ir  ii  to 
<VM*totototototo 
I  •••••• 

•  4  W  I  I  I  I  I  I 


O  PI  to  C  *4  2  4 


ON J «CUN 


conoono 
u  o  n  n  o  o  b 
t>  o  o  o  o  o  e 

OOO0OOO 


«UMWWWUW 


tr 

s 

tC 


u**»»«*« 

WNNNNNNN 


0O  4  •  ♦  •  *  «  • 

to  m  Ui  ui  ui  ui  ut  UJ  ut 

ut  to  •  tr  ft-  to  cr  <r  *r 


«  UJ  I  I  M  I  •  • 


o  m  N  i  to  n  « 
ON7JOO  N 


to  to  to  to  to  to  to 


7-37 


Fig.  7-15  (Cont.) 
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Fig.  7-15  (Cont.) 
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7.6  SAMPLE  CASE  6  -  ELLIPSOID 


The  critical  external  pressure  according  to  bifurcation  theory  is  to  be  determined  for 
an  orthotropic  ellipsoid  shell  with  clamped  circular  boundaries.  The  geometry  of  the 
shell  is  shown  in  Fig.  7-16. 

Analysis  with  BOSOR4  (Ref.  15)  shows  that  for  this  shell  the  critical  external  pressure 
is  4327  psi  for  a  mode  with  2  circumferential  waves.  Therefore,  a  shell  covering 
1/8  of  the  circumference  can  be  analyzed. 

The  meridional  boundary  conditions  are  the  same  as  for  Case  4.  The  input  cards 
associated  with  this  case  are  displayed  in  Fig.  7-17.  Portions  of  the  output  are 
presented  in  Fig.  7-18. 

The  Critical  external  pressure  determined  here  is  4741  psi.  The  discrepancy  in 
critical  external  pressure  is  due  to  the  coarseness  of  the  grid  used  here. 
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Fig.  7-16  Sample  Case  6  -  Orthotropic  Ellipsoid 


SAMPLE  CASE  6 
INPUT 


SAMPLE 

CASE  6  -  ELLIPSOID 

ORTHOTROPIC 
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2  20 
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Fig.  7-17  Display  of  Input  Cards  for  Sample  Case  6 
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SAMPLE  CASE  6  -  OUTPUT 
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Fig.  7-18  Excerpt  of  Output  for  Sample  Case  6 
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7.7  SAMPLE  CASE  7  -  HYPERBOLOID 


The  critical  axial  line  load  according  to  bifurcation  theory  is  to  be  determined  tor  a 
hyperboloid  shell  segment  clamped  at  the  large  circular  edge  and  constraint  in  all  but 
the  axial  direction  at  the  other  edge.  The  geometry  of  the  shell  is  shown  in  Fig.  7-22. 

Analysis  with  BOSOR4  (Ref.  IS)  shows  that  for  this  shell  the  critical  axial  line  load 
is  1740  lb/ in.  for  a  mode  with  6  circumferential  waves.  Therefore,  a  shell  covering 
1/24  of  the  circumference  can  be  analyzed. 

The  meridional  boundary  conditions  are  the  same  as  for  Sample  Case  4.  The  input 
cards  associated  with  this  case  are  displayed  in  Fig.  7-23.  Portions  of  the  output  are 
presented  in  Fig.  7-24. 

The  critical  axial  line  load  determined  here  is  1927  lb/in.  The  discrepancy  in  the 
critical  load  is  due  to  the  coarseness  of  the  grid  in  the  longitudinal  direction,  selected 
here  to  reduce  the  run  time. 
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Axis  Of 
Symmetry 


The  Equation  Of  The 
Generator  is: 


And  R*R,  ■  1.0 
0*  15* 


Fig.  7-19  Sample  Case  7  -  Hyperboloid 
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Fig.  7-20  Display  of  Input  Cards  for  Sample  Case  7 
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Fig.  7-21  Excerpt  of  Output  Sample  Core  7 
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Fig.  7-21  (Cont.) 


CO  3  C  •O.tflM'C-Ol  -5.8S05O0C-01  (.  •l.tlMtiMI  5.05777»C-0i  S. 

•MI  -MtlMU-lt  *5*  0$  054  OE-0 1  -MJJlfPM*  -MltlMMt  5. 057  ?7<.C- 0»  (.07(177 

MOO  -5.0  5  I  78  9C*0 1  -5.00051.01-01  -1.N13751C-10  - 1.8  1  7<.7  0C-0(  5. 057 77SC- 04  I.I'.OS'IO 


7-53 


THf  (UM1.IN5  1.000  0*570  ON  l  INC  OK  OICVUCtTION  7HC0NY  IS  MtlMKHl  TINES  TNt  STOKTINS  1000 


a 

w 

v 


.  K  N  K  K  ( 


K  N  K  H  »)  n  r  K 

MNN4  ~  — 

•  I  I  •«  I  I  I 


y  W  wwwiuu^ 
►  C,n*i'^e\o/ 
«c(r-/ooi/'('j\ 

CN"O»0«»K 

^  J  himt  o**o 


uuyuw u  k  u 

►  fs/^«irc“«=»oiA 

hlAV'tfClT  *<  »  H 
C*  •»  *k  k.  k>  *«  *0  .4  O 

oo«*-<r^»oo^  .* 

^KlNHff-erjk. 


KKKM>.K 

I  I  I  •  •  • 

WkiWUWW 
tft  Mr  Aj  N  M 
N  4  N  <0  *0  <0 
r  MT  <U  N  N 
«  N  J  N  N  N 
<300000 
N  N  ifi  if  X 

I  I  I  •  •  • 


I  I  I  I  I  <  I  I 
UUWUWWW w 
K  nmK40tf)0<0 

X  ►■|f>0OO0NW 

k-  r^<0»w 

hi  BtThOO^  HK 

o  4rjoN»(»«i0 


•  •»•••#• 


OOOOOOOO 

I  I  I  I  I  I  •  • 
yuuuuwuw 
n  IT  K  #  o  «  « 

3  MT^OO^N- 
Khoir»r^N 
cr^ori  •<  N 


I  I  •  •  I  I  •  • 


OOOOOOOO 

t  I  I  I  I  I  I  I 
hlhUhiibWuy 
H 

k-  OK.fNitf>—  *0  *\J  <0 

c  ^Ko-iroiir 


OOOOOOOO 


i  •  • 


i  I 


•  ♦♦♦*♦♦♦ 
WWUWWUUU 
x  00inr,i^4N 

k-  trs<r«\»«:rv*ki  —  *0 

Ui  «*»(V<C*^*»oa,irtO 

< mw4*>*+^**T<r 


*£>&'*••***•  —  XO 
U'lWJK'ffCK 

j  j  ir**  j  ir-c 
ecooiror«0 

*C>  O  K  ®  O' 


OOOOOOOO 

I  t  I  •  I  »  I  • 

wuuwuuuy 
i r  0  «fl  o  n  mt 
oiriH^^ow 
k>  k.  k.  0  O  >1  C  K 
O0k0«HlTK 
O  4  0  N  iflrrK 


OH0M000KN 


I  I  I  •  <  • 

ItlWkiWUlti 
I  O  O  O  0>  0  o 
r  cmto 

khh  4  40 

t  tr  tr  O'  O'  O'  *. 

>  0  0  0  V)  10 

00«nno 

OMNN00HOD 

•  I  •  .1  •  • 

i 

K  K  K  K  K  k.  3 


•  ••••• 

towbwww 

3  k  k  k  IT  1  »< 
k  k  k  0  ^  H 
k  k  k  0  0  H 

•  00NK0 
0»  m  wt  k  x  tf\ 

OHHMNN0C 

I  f  I  •  I  • 


•  I  I  I  «  •  I  • 

0000004^1 


OOOOOOOO 
I  I  f  I  I  f  I  I 
WU«UI»1  w  Uiwb 
O0  (Tk  Hiprtk 

te*<o000oc 

OkP0OO-fj 
ONXO404O 
M  r  ►•  «  o  r  ^  k 
^r>n«4O04k 


0H0P10kkN 

>^0^040^0 

040^0000 

OH0N0 00k 
0  0  N  «  N  0  k  0 


Mkivwwi, 


i  o  o  o  o  o  o 


OOOOOOOO 
<  I  I  «  I  I  I  • 
UWWb WWbU 
4N0400N0 
:  o^»-44'4'o»^o 

4  N  N  4  *■  »1 
«■  v  k  n  k  .<u  u 
0N0k.«4-*O 
«T  OOOOIMCO 

O  »4  n  0  k  *  0  0  »4 


J  UJ  U  U;  U>  U  W  U 
.  y  -r  CTK  *lk 
^  0  o  4  0  0  0  0 
<  k  >i  ”  4  o  if  a 
i  *■*  k.  <0  w'  *n  -*0 
i  4  i  C  .*>.^1  e 
r  O  4»  CT  o  *0  k.  k. 


i  ^NNNf 


0000000004 


>0O0O0O0rt 
»  k  0  rj  o  k  4  k  n 
)  0  k  0  0  n  kl  ••  O 


>  0  n  o  o  a  r  o  o 

'000000  4  0 

>  N  O  k  0  f  j  O 

•  <0  ii>  r»  k.  —  O 


ooooooono 
000000000 
.  fiki/'kr)k4kO 
(JC^>00Hf|xO 


O  0  O  0  O  0  O  | 

:  ►  o  * 
n  «  » 


7-54 


Fig.  7-21  (Cont.) 


7.8  SAMPLE  CASE  8  -  PARABOLOID 


The  critical  external  pressure  according  to  bifurcation  theory  is  to  be  determined  for 
a  paraboloid  shell  with  clamped  circular  edges.  The  geometry  of  the  shell  is  shown 
in  Fig.  7-22. 

Analysis  with  BOSOR4  (Ref.  15)  shows  that  for  this  shell,  the  critical  external  pressure 
is  714  psi  for  a  mode  with  10  circumferential  waves.  Therefore,  a  shell  covering  1/40 
of  the  circumference  can  be  analyzed. 

The  meridional  boundary  conditions  are  the  same  as  for  Sample  Case  4.  The  input 
cards  associated  with  this  case  are  displayed  in  Fig.  7-23.  Portions  of  die  output  are 
presented  in  Fig.  7-24. 

The  critical  external  pressure  determined  here  is  703  psi. 
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pig.  7-24  Excerpt  of  Output  for  Sample  Case  8 
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7.9  SAMPLE  CASE  9  -  FIBERWOUND  CYLINDER 

The  critical  axial  line  load  according  to  bifurcation  theory  is  to  be  determined  for  a 
circular  cylindrical  shell  built  up  from  3  fiber-wound  layers  and  clamped  at  both  ends. 
The  geometry  of  the  shell  is  shown  in  Fig.  7-25. 

Analysis  with  BOSOR4  (Ref.  15)  shows  that  for  this  shell,  the  critical  axial  line  load  is 
500  lb/in.  for  a  mode  with  11  circumferential  waves.  Therefore,  a  shell  covering  one 
half  the  length  and  1/44  of  the  circumference  can  be  analyzed. 

The  meridional  boundary  conditions  are  the  same  as  for  Sample  Case  4.  The  input 
cards  associated  with  this  case  are  displayed  in  Fig.  7-26.  Portions  of  the  output  are 
presented  in  Fig.  7  -27. 

The  critical  axial  line  load  determined  here  is  609  lb/in.  Note  that  the  discrepancy  in 
the  solution  is  due  to  the  coarseness  of  the  grid  in  the  longitudinal  direction,  selected 
here  to  reduce  the  run  time. 
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Fig.  7-25  Sample  Case  9  -  Fiberwound  Cylinder 
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Fig.  7-27  Excerpt  of  Output  for  Sample  Case  9 
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7.  10  SAMPLE  CASE  10  -  THERMAL  LOADS 


The  critical  load  combination  according  to  bifurcation  theory  will  be  determined  for 
a  simply  supported  circular  cylindrical  shell  loaded  with  axial  compressive  line  load 
and  subjected  to  temperature  distribution  that  varies  linearly  along  the  length  and  is 
constant  circumferentially. 

Analysis  with  BOSOR4  (Re  15)  shows  that  for  this  shell,  the  lowest  eigenvalue  of 
3. 156  occurs  in  conjunction  with  26  circumferential  waves.  Therefore,  a  shell  covering 
one  half  the  length  and  1/52  of  the  circumference  can  be  analyzed. 

User-written  subroutine  MATER  is  used  here  to  introduce  the  material  properties 
and  the  temperature  distribution.  However,  this  particular  example  could  be  solved 
using  user-written  subroutine  TEMP  and  input  card  M-2B.  The  geometry  and  tempera¬ 
ture  variation  of  the  shell  are  shown  in  Fig.  7-28.  Subroutine  MATER  and  the  input 
cards  associated  with  this  case  are  shown  in  Fig.  7-29  and  Fig.  7-30,  respectively. 
Portions  of  the  output  are  presented  in  Fig.  7-31.  The  critical  load  for  this  case  is 
the  eigenvalue  times  (STLD*PP),  where  PP  is  the'base  load  composed  of  the  axial 
line  load  (1  lb/in. )  plus  the  thermal  loads  due  to  the  temperature  distribution  as  shown 
in  Fig.  7-28.  The  lowest  eigenvalue  determined  here  is  3. 269. 


AT  =  100 


Fig.  7-28  Sample  Case  10  -  Thermal  Loading 
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:  jh  '  »*>*<» 


rU!rOUTIf;€  MATER  <  X  ,  Y  ,  IP  ,  TO  E  G ,  £  X  ,  E  Y  ,  U,G  »  A 1 ,  A  ?) 

DIMENSION  TDEG<  IP)  ,EX(IP),EY(IP),U{IP),G(IP>,A1(IP>,A2(IP> 

rnr^ON  /  OFS  T  /  TC.Z 

c  t o  rn  i  must  me  set  by  the  user 

C  1C  =  TOT/L  THICKNESS  OF  SHELL  (TO  =  AT) 

C  7  -  DISTANCE  FROM  REFERENCE  SURFACE  TO  MIDSURFACE  OF  SHELL  WALL 

A  T  =  :  .03 

z  =  o . 
n-  at 

ro  l  l  =  l  *  ip 
c--i . 

TCEfCL)  =  6 0.*  (2.-X/10.) 
rx(D  =c*i  ooo  noon. 

fY{L)=C#10000GQ0. 

L!  ( L  )  =  3 .  3 
G  (L  )=C*4CQ0  0  00. 

AI(  L) =0 . COO  1 
A2 ( L )  =  1 •  0  00  1 
1  CONTINUE 
RETURN 
FND 


Fig.  7-29  Subroutine  MATER  for  Sample  Case  10 


SAMPLE  CASE  10 
INPUT 


SAMPLE  CASE  1C  -  THERMAL  LOADS  (IWALL=8> 
ill 
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G-2 
D-l 
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B-l 
L-l 
L-2 
P-1  A 
0-1 
0-2 
M-l 


Fig.  7-30  Display  of  Input  Cards  for  Sample  Case  10 
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SAMPLE  CASE  10  -  OUTPUT 
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7-31  Excerpt  of  Output  for  Sample  Case  10 
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Fig.  7-31  (Cont.) 


4 


7.11  SAMPLE  CASE  11  -  PLASTICITY 

The  stress  distribution  in  the  interior  of  a  square  plate  subjected  to  concentrated  shear 
force  at  one  edge  is  to  be  determined  using  the  plasticity  branch  of  the  STAGS  computer 
program.  The  geometry  and  boundary  conditions  are  shown  in  Fig.  7-32. 

The  stress-strain  curve  is  described  in  the  input  data  utilizing  five  distinct  points 
(material  components).  The  equivalent  White-Bcsscling  material  components  are 
shown  in  the  output.  Since  the  stresses  are  uniform  across  the  thickness  in  this 
problem,  only  three  points  are  used  for  the  numerical  integration  through  the  thickness 
of  the  plate. 

Points  at  which  the  effective  stress  is  in  the  inelastic  range  are  designated  with  an 
asterisk  as  shown  in  the  output.  To  continue  a  nonlinear  plasticity  nm,  only  the  last 
record  can  be  used  for  a  restart  (ISTART=3),  because  the  strains  for  records  1  and  2 
are  not  saved . 

The  input  cards  associated  with  this  case  are  displayed  in  Fig.  7-33,  Portions  of  the 
output  are  presented  in  Fig.  7-34. 
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Concentrated  Shear 
Force  at  Mid  Length 


Fig.  7-32  Sample  Case  11  -  Plasticity  (Plate) 
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Fig.  7-33  Display  of  Input  Card  for  Sample  Case  11 
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Fig.  7-34  Excerpt  of  Output  for  .p  e  11 
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Section  8 

STAGS  POST  PROCESSOR 


8.1  INTRODUCTION 

The  STAGS  Post  Processor  (STAGPP)  is  designed  to  be  used  with  the  STAGS 
program.  The  function  of  the  STAGS  Post  Processor  is  to  present  the  cal¬ 
culated  and  stored  data  of  STAGS  in  a  two  dimensional  or  contour  plot  out¬ 
put.  The  STAGS  Post  Processor  tekes  the  data  stored  on  tape  19  by  the  STAGS 
program  and  creates  a  printer  plot,  a  CAL-COMP  pen  plot,  or  both  depending 
on  the  requirements  of  the  iser.  The  STAGS  Post  Processor  is  designed  to 
plot  selected  output  and/or  contour  plots  of  displacements  (w,  v,  u,  p^,  p^), 
and  stress  resultants  (Nx>  N^,  M^,  M^,  M^)  .  Also,  the  load  displace¬ 

ment  history  and  collapse  mode  plots  are  included  in  the  post  processor's 
capabilities . 

This  version  of  the  STAGS  Post  Processor  was  written  for  the  CDC  6600  computer 
utilizing  a  model  763,  zip  mode  10"  drum  CAL-COMP  pen  plotter.  The  user  must 
supply  the  appropriate  computer  system  controls  when  a  CAL-COMP  plot  is 
requested . 

There  are  two  main  branches  of  the  STAGS  Post  Processor.  The  first  branch 
allows  for  two  dimensional  plots  of  displacements,  stress  resultants  and 
load  displacement  history  and  the  second  branch  does  only  the  contour  plot6 
of  displacements,  stress  resul+^t.s,  and  the  collapse  mode.  Figure  8-1  shows 
a  typical  data  deck  format.  Table  8-1  is  a  summary  of  the  input  cards. 
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Linear  Solution 


Selected  Displacements  with  Editing 


Selected  Stress  Resultants  with 
Editing 


Load  Displacement  History  - 
Nonlinear 

Collapse  Mode  Contour  -  Nonlinear 


Special  Output  - 
Displacement  Contour 


Special  Output  - 
Stress  Resultant  Contour 


Special  Output  - 
Stress  Resultant  Edited 


Normal  Terminator  for  STAGS 
Post  Processor 


Fig.  8-1  STAGS  Post  Processor  Sample  Data  Input  (Sample  Case  3>  Page  7-22) 


8-2 


Table  8-1 

STAGS  POST  PROCESSOR 


SUMMARY  OF  INPUT  CARDS 

Card  Symbol  Format 

PI  IPLT,  IPLTSW  l6l‘o 

Include  Card  P2  only  if  0  £  IPLT  s  3  or  IrLT  -  '(  or  8.  If 
^  i  IPLT  s  6  go  to  the  P7  card. 

P2  I  AXIAL,  I  OUT,  IEDIT,  IHST  l6l  5 

Include  cards  P3  through  P6  only  if  IEDIT  =  1  on  the  P2  card. 

P3  NSTEP,  NXYEDT  l6l'; 

PU  JSTEP(I),  1=1,  NSTEP  1615 

P5  XYEDIT(I),  1  =  1,  NXYEDT  8F10.S 

Pb  IVAR(I),  1=1,  NXYEDT  8F10.5 

Pb  IVAR(I),  1=1,6  l6lb 

Include  cards  P7  through  PIO  only  if  ^  s  IPLT  s  6  on  the  PI  card. 

P7  NROW ,  NCOL,  NPB ,  NCTOUT,  NRW1,  NRW2,  NCL1  1615 

P8  NB(l),  1=1,  NPB  1615 

P9  IVAR(l),  1=1,6  l6l 5 

PIO  HEAD(I),  1=1,3  3A10 


Include  as  many  PIO  cards  as  there  are  IVAR(l)  =  1  on  the 
P9  card. 


8-3 


8.2  INPUT  DESCRIPTION 


PI  Plot  Data  Card 

This  card  is  used  to  select  how  and  which  stored  data  is  to  be  plotted. 


Variable 

Format 

Columns 

Description 

IPLT 

15 

1-5 

IPLT  selects  from  a  file  created  by  STAGS 
the  data  which  is  to  be  plotted. 

IPLT=0  Plot  linear  solution. 

IPLT=1  Plot  selected  displacements. 

IPLT=2  Plot  selected  stress  resultants. 

IPLT=3  Plot  load  displacement  history 
(nonlinear  analysis  only) 

IPLT=4  Plot  collapse  mode  (nonlinear  only) 

IPLT=5  Displacements  contour  plot 

(special  output  see  Page  6-33) 

IPLT =6  Stress  resultants  contour  plot 
(special  output). 

IPLT=7  Displacements  plot  (special  output). 

IPLT=8  Stress  resultants  plot  (special 
output ) . 

IPLT=999  Normal  terminator  for  STAGS  Post 
Processor . 

IPLTSW  15  6-10  IPLTSW=1  Do  both  printer  and  CAL-COMP 

pen  plots. 

IPLTSW=2  Do  only  printer  plots. 

IPLTSW=3  Do  only  CAL-COMP  pen  plots. 

IPLTSW=999  Normal  terminator  for  STAGS 
Post  Processor. 

NOTE:  If  It  s  IPLT  s  6  go  to  the  P7  card. 


8-4 


P2  Select  Data  Card 


i 


< 

r 

t 

j 

i 

The  select 

data  card 

allows  the  user 

to  exercise 

certain  options  concerning 

j 

which  load 

to  select, 

print  output, 

editing  capability  and  load  displacement 

history. 

j 

! 

i 

Variable 

Format 

Columns 

Description 

{ 

I AXIAL 

15 

1-5 

IAXIAL=0 

Plot  load  factor,  PA. 

0 

IAXIAL=1 

Plot  load  factor,  FB. 

IOUT 

15 

6-10 

I017T=0 

Do  not  print  plotted  data. 

i 

I017T=1 

Print  plotted  data. 

i 

* 

IEDIT 

15 

11-15 

IEDIT =0 

Plot  everything  requested  on 
the  PI  card.  No  more  cards 
are  needed. 

i 

IEDIT =1 

Plot  only  selected  or  special 
output  requested  on  the  P3  to 

P6  cards. 

IHST 

15 

16-20 

IHST=0 

Do  not  plot  a  load  displacement 
history . 

IHST=1 

If  IPI/T=3  do  a  CAL -COMP  load 
displacement  history  pen  plot. 

P^^di^Control  Card 

Used  only  if  IEDIT*1  on  the  P2  card.  The  edit  control  card  states  how  many 
load  steps  and/or  rows  or  columns  to  select  for  plotting. 


Variable 

Format 

Columns 

Description 

NSTEP 

15 

1-5 

NSTEP  is  the  number  of  load  steps  requested. 
If  no  load  steps  are  requested  set  NSTEP*1 
and  see  card  P^ . 

NXYEDT 

15 

6-10 

NXYEDT  is  the  number  of  rows  or  columns 
requested.  If  no  rows  or  columns  are 

requested  set  NXY1DT=1  and  see  card  P5. 

| 

t 
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Used  only  if  IEDIT«1  on  the  P2  card, 
step  requested  for  plotting. 

f 

> 


l 

Variable 

Format 

Columns 

i 

* 

* 

JSTEP(I) 

1=1,  NSTEP 

1615 

l-8o 

l 

t- 

> 

: 

s 


This  card  specifies  the  individual  load 


Description 

Actual  load  step. 

JSTEP(l)=777  No  editing  by  load  step 
performed  (use  for  linear 
and  bifurcation  analysis). 

JSTEP(l)=  Actual  load  step  requested 
(use  for  nonlinear  analysis 
only) . 


P^^Selec^Roj^^Co^am^ard 

Used  only  if  IEDIT=1  on  the  P2  card.  This  card  specifies  the  individual  rows  and 
columns  requested  for  plotting. 


Variable 

Format 

Columns 

Description 

XYEDIT(l) 

8F10.5 

i-8o 

Actual  rows  or  columns . 

I=1,NXYEDT 

XYEDIT(l)  =777-0  No  editing  by  rows  or 

columns  performed. 

XYEDIT(l)  =  Actual  coordinate  of  rows 

or  columns  requested. 

) 

t 

i 

j 

P6  Select  Displacement  or  Stress  Resultant  Card 


Used  only  if  IEDIT=1  on  the  P2  card.  This  card  reads  data  into  an  array  which 
indicates  which  displacements  or  stress  resultants  are  required.  Displacements 
are  requested  if  IP1/T=1  or  7  on  the  PI  card.  Stress  resultants  are  requested 
if  IFLT=2  or  8  on  the  PI  card. 
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Variable 

Format  Columns 

Description 

IVAR(I) 

1=1,6 

615  1-30 

IVAR(I)=0  Do  not  plot 

IVAR(l)=l  Plot 

1=1 

W  displacement  or  NX  stress 

resultant 

1=2 

V  displacement  or  NY  stress 

resultant 

1=3 

11  displacement  or  NXY  stress 

resultant 

1=4 

BX  displacement  or  MX  stress 

resultant 

1=5 

BY  displacement  or  MY  stress 

resultant 

1=6 

MXY  stress  resultant 

NOTE:  Include  cards  P7  through  P10 

only  if 

4  £  IPLT  s  6  on  the  PI  card. 

P7  Contour 

Control  Card 

The  contour  control  card  allovs  the  user  to  define  the  grid  for  contour  plots. 
NROW  and  NCOLS  (NR  and  NC  in  ST^GS)  define  the  grid  and  NPB  is  the  number  of 
points  defining  the  plot  boundaries  stored  in  the  array  NB.  The  NROW  and 
NCOL  define  the  NB  grid  points  vhere  1  si  s  NROW  and  1  s  J  s  NCOL.  Then 
NB  =  (l-l)*NCOL)  +  J  (see  Fig.  8-2).  If  a  cutout  is  present  NB  will  not 
include  the  grid  points  of  the  cutout  (see  Fig.  8-3).  Once  the  NB  array 
has  been  established  the  user  can  select  a  section  of  the  total  shell  by 
choosing  the  proper  NB  values.  Fig.  8-2  where  (NB(l),  1*1,  NPB)  =1,9, 

14,  12,  7  will  create  a  contour  plot  of  the  center  portion  of  the  flat  plate 
only. 

If  IPLT=6  on  the  PI  card  and  NCOL  >  NROW  on  the  P7  card,  refer  to  Figs .  8-4 
and  8-?  as  examples  of  NCTOITT,  NRW1,  NRW2,  and  NCL1. 
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Variable 

Format 

Columns 

Description 

NROW 

15 

1-5 

Number  of  rows  (stations  along  the  X  coordinate). 

NCOL 

15 

6-10 

Number  of  columns  (stations  along  the  Y 
coordinate ) . 

NPB 

15 

11-15 

Number  of  points  stored  in  array  NB  which 
define  the  plot  boundary. 

NCTOUT 

15 

16-20 

Number  of  grid  points  of  the  cutout  to  be 
excluded  (see  Fig.  8-3). 

NRW1 

15 

21-25 

Row  number  of  one  edge  of  cutout  (same  as  in 
the  STAGS  program) . 

NRW2 

15 

26-30 

Row  number  of  other  edge  of  cutout  (same  as 
in  the  STAGS  program) . 

NCL1 

15 

31-35 

Column  number  of  edge  of  cutout  (same  as  in 
the  STAGS  program) 

P8  Plot  Grid  Card 

This  card  reads  the  MB  numbers  which  specify  the  contour  plot  boundary. 
Successive  pairs  of  points  listed  in  NB  will  be  Joined  by  a  straight  line. 
Note  the  first  point  in  NB  should  be  repeated  as  the  last  point  in  NB  in 
order  to  close  the  contour  plot. 


Variable 

Format 

Column 

D  ocriptior. 

NB(I) 

1=1,  NPB 

1015 

■*  ;o 

The  actual  NB  points  which  define  the 
plot  boundary  (see  Figs.  8-2  through 
8-5). 

Select  Displacement  or  Stress  Resultant  Card 


This  card  reads  data  into  an  array  which  indicates  which  displacements  or 
stress  resultants  are  required.  If  IPLT=5  on  the  PI  card,  we  have  a  con¬ 
tour  plot  of  displacement.  If  IPLT=6  on  the  PI  card,  we  have  a  contour 
plot  of  stress  resultan4  . 
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Variable 

Format 

Columns 

Description 

IVAR(l) 

615 

1-30 

IVAR(I)=0 

Do  not  plot 

1=1,6 

IVAR(I)=1 

Plot 

1=1  W  displacement  or  HX  stress  resultant. 

1=2  V  displacement  or  NY  stress  resultant. 

1=3  U  displacement  or  NXY  stress  resultant. 
1=4  MX  stress  resultant. 

1=5  MY  stress  resultant. 

1=6  MXY  stress  resultant. 


PIO  Plot  Heading  Card 


The  user  must  read  individual  titles  for  each  contour  plot.  There  are  8s 
many  PIO  cards  as  there  are  IVAR(l)  requested  on  the  P9  card. 


Variable 

Format 

Colums 

Description 

HEAD (I ) 

3A10 

1-50 

User  supplied  heading  for  contour  plot 

1=1,3 

8-9 


ROW  1 


ROW  2 


ROW  3 


ROW  4 


COL  1  COL  2  COL  3  COL  4  COL  5 
12  3^5 


6 

7 

8 

9 

11 

14 

16 

17 

18 

19 

10 


15 


20 


NR  <  NC 


NR0W=4,  NC0L=5,  NPB=5,  NCTOUT=0,  NRW1=0,  NRW2=0,  NCL1=0,  (NBClJ^l.NPB)  = 
[NB:1, 5,20,16,1 } 

Fig.  8-2  Flat  Plate  Contour  Control  Card  Setup 


COL  1  COL  2  COL  3  COL  4  COL  5 


ROW  1 

1 

2 

3 

4 

5 

_  Y 

ROW  2 

6 

7 

8 

9 

10 

NR  <  NC 

ROW  3 

11 

12 

13 

14 

15 

ROW  4 

^  - 

16 

17 

18 

19 

I 

X 


NROW-4,  NC0L=5,  NPB=7,  NCT0UT=1,  NRW1=3,  NRW2=4,  NCL1=2,  (NB(l),I=l,NPB)  = 
{NB:1, 5,19,16, 12, 11,3  } 

Fig.  8-3  Flat  Plate  with  Cutout  Contour  Control  Card  Setup 
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(STAGS  Poet 
Processor) 


(STAGS) 


COL  1 
ROW  1 


COL  2 
ROW  2 


COL  3 
ROW  3 


COL  4 
ROW  4 


(STAGS  Post  Processor) 
(STAGS) 


ROW  1 

COL  1 

1 

2 

3 

ROW  2 

COL  2 

5 

6 

7 

ROW  3 

COL  3 

9 

10 

11 

ROW  4 

COL  4 

13 

14 

15 

ROW  5 

COL  5 

17 

18 

19 

NR0W=4, 

NC0L=5,  NPB=5 

Y 

,  NCTOUT 

=0,  NRW1  = 

0,  NRW2=0 

16 


20 


NR  <  NC 


{NB:1,U, 20,17,1} 


Fig.  8-4  Flat  Plate  Contour  Control  Card  Setup  When 

NCOL  Greater  Than  NROW  (Stress  Resultants  Only) 


(STAGS  Post 

UUL  1 

COL  2 

COL 

3 

COL 

Processor) 

(STAGS) 

ROW  1 

ROW  2 

ROW 

3 

ROW 

ROW  1 

COL  1 

1 

2 

3 

ROW  2 

COL  2 

4 

5 

6 

1  " 

1 

1 

i  7 

ROW  3 

COL  3 

8 

9 

10 

ii 

ROW  4 

COL  4 

12 

13 

14 

15 

ROW  4 

COL  5 

16 

17 

18 

19 

Y 

(STAGS  Post  Processor) 
(STAGS) 

X 


NR  <  NC 


NOTE:  NRW1,  NRW2,  NCL1  are  in  the  STAGS  Post  Processor  Notation 

NR0W=4,  NC0L=5,  NPB=7,  NCT017T  =  1,  NRW1=1.  NRW2=2,  NCL1=3,  (NB(l).I  1,(NPB)  -= 

[MB:1, 3, 6, 7,19,16,11 

Fig.  8-5  Flat  Plate  \ith  Cutout  Contour  Control  Card  Setup 

Wher.  NCOL  Greater  Thar.  NROW  (Stress  Resultants  Only) 
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